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ABSTRACT 


An  optimization  analysis  is  presented  for  nozzles  with  gas-particle 
flows.  The  problem  is  formulated  to  maximize  the  axial  thrust  produced  • 
along  the  nozzle  contour  for  a  general  isoperimetric  constraint  such  as 
constant  nozzle  length  or  constant  nozzle  surface  area.  The  effects  of 
the  ambient  pressure  are  included  in  the  thrust  expression  to  be  maximized. 
The  characteristic  and  compatibility  equations  are  developed  and  numerical 
techniques  are  presented  for  use  in  conjunction  with  the  characteristic 
and  compatibility  equations.  A  solution  procedure  is  presented  which 
determines  whether  or  not  a  given  nozzle  contour  is  an  optimal  solution 
and  a  relaxation  technique  is  presented  which  adjusts  the  nozzle  contour 
toward  the  optimal  solution.  Selected  parametric  studies  are  presented. 
These  studies  illustrate  the  effects  of  changing  mesh  size,  particle 
size,  particle  mass  flow  rate,  inlet  angle,  drag  coefficients,  heat 
transfer  coefficients,  throat  radius  of  curvature,  and  the  scale  on  the 
thrust  performance  and  the  nozzle  geometry  of  the  optimal,  fixed  length 
nozzle. 


FOREWORD 


The  present  study  is  part  of  the  program  "An  Analytical  Study  of 
the  Exhaust  Expansion  System  {Scramjet  Scientific  Technology)'*  being 
conducted  by  the  Jet  Propulsion  Center,  Purdue  University,  Lafayette, 
Indiana,  under  United  States  Air  Force  Contract  No.  F33615-67-C-1068, 
Project  3012,  Task  301209,  BPSN  7(63  301206  62G5214).  The  Air  Force 
program  monitor  was  Capt.  Gary  0.  Jungwirth  of  the  Air  Force  Aero 
Propulsion  Laboratory  (AFAPL/RJT).  This  report  presents  the  formula¬ 
tion,  numerical  solution  procedure  and  the  results  of  selected  para¬ 
metric  studies  of  the  design  of  maximum  thrust  nozzles  with  gas- 
particle  flows.  Volume  II  is  the  computer  program  user's  manual. 

Tnis  report  was  submitted  by  the  authors  on  31  May  1971. 

Publication  of  this  report  does  not  constitute  Air  Force  approval 
of  the  report's  findings  or  conclusions.  It  is  published  only  for  the 
exchange  and  stimulation  of  ideas. 
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SECTION  I 


INTRODUCTION 

Many  propellant  confcinations  produce  condensed  phases  In  the  exhaust 
products.  These  condensed  phases  may  be  due  to  the  introduction  of  metal 
additives  designed  to  increase  energy  release.  These  condensed  phases, 
however,  introduce  performance  losses  due  to  the  non-equilibriirs  effects 
of  heat  transfer  and  drag  between  the  gas  and  the  condensed  phase 
particles. 

The  first  application  of  optimization  techniaues  to  the  design  of 
rocket  nozzles  v?as  made  by  Gudeney  and  Hanisch  (1)  for  hosaentropic  flew 
in  1955.  Rao  (2,3)  simplified  tfc«  analysis  and  applied  the  formulation 
of  the  problem  to  standard  nozzles  and  to  plug  nozzles-  Suderiey  (4) 
then  extended  the  results  to  isentropic  flows  which  allow  entropy  to 
vary  between  streamlines. 

Figure  1  reoresents  the  general  sadel  used  for  formulating  the  optim¬ 
ization  of  standard  axi symmetric  nozzles.  Iri  the  above  analyses,  the 
problem  was  fbmilated  to  provide  maximum  thrust  across  an  exit  control 
surface,  BC,  and  employed  a  constant  length  desigi  constraint.  Since  the 
formulation  is  limited  to  the  exit  control  surface,  no  dissipative  effects 
in  the  flaw  field  are  allowed  and  nozzle  design  constraints  not  directly 
related  to  the  exit  control  surface  are  not  available. 
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FIGURE  1.  GENERAL  OPTIMIZATION  MODEL 


Guderley  and  Armitage  (5,6)  formulated  the  optimization  problem  over 
the  entire  region  ABC  in  order  to  use  a  constant  surface  area  as  a  design 
constraint  in  lieu  of  constant  length.  This  approach  permits  the  use  of 
a  wide  range  of  geometric  design  constraints.  The  complexity  of  the 
formulation  and  the  numerical  solution  are  greatly  increased  over  the  pre¬ 
vious  formulations. 

The  Guderley-Armitage  approach  can  also  be  extended  to  dissipative 
flows  since  the  entire  region  ABC  is  considered  in  the  formulation. 

Hoffman  and  Thompson  (7,8)  formulated  the  problem  for  gas-particle  flows 
and  Hoffman  (9,10)  formulated  the  problem  for  reacting  non-equilibrium 
flows.  Further  work  was  performed  at  Purdue  University  to  develop  the 
numerical  schemes  and  to  furnish  working  computer  programs  for  the  design 
of  maximum  thrust  nozzles  having  flow  fields  with  dissipative  effects. 
Scofield  and  Hoffman  (11)  treated  rotational  or  non-equilibrium  simple 
dissociating  gas  flows,  Humphreys,  Thompson  and  Hoffman  (12)  treated  plug 
nozzles  with  fixed  inlet  geometry,  and  Johnson,  Thompson  and  Hoffman  (13) 
treated  plug  nozzles  with  variable  inlet  geometry. 

This  work  presents  the  formulation  and  nisnerical  schemes  developed 
to  determine  thrust  r.czzle  contours  for  nozzles  with  condensed 

particles  in  the  flow  field,  and  presents  the  results  of  an  extensive 
parametric  study.  The  optimization  problem  is  formulated  over  the  region 
ABC  (Fig.  1),  and  follows  that  presented  in  Refs.  (7,8).  The  method  pre¬ 
sented  herein  uses  the  calculus  of  variations  to  develop  Lagrange  multi¬ 
plier  equations,  develops  a  numerical  scheme  to  apply  these  equations  to 
a  previously  calculated  flow  field  in  an  assumed  nozzle  contour,  deter¬ 
mines  an  error  function  along  the  assumed  wall,  and  calculates  a  new  wall. 
This  procedure  is  repeated  until  the  error  function  goes  to  zero,  and  the 


rtisvimum  thrust  rnntot’i'  is  obtained.  An  existing  gas-particle  flow  field 
analysis  program  developed  by  Kliegel  and  Nickerson  (14)  was  adapted  to 
provide  the  evaluation  of  the  flow  properties  in  the  nozzle. 
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SECTION  II 
ANALYSIS 


1 .  INTRODUCTION 

In  formulating  thrust  optimization  problems,  two  basic  approaches 
have  been  used.  The  first,  and  by  far  the  easiest  to  use,  maximizes  the 
thrust  written  in  terms  of  flow  variables  axross  an  exit  control  surface. 
However,  the  application  of  this  approach  is  limited  since  the  formula¬ 
tion  along  the  exit  surface  does  not  allow  for  dissipative  flows  or  for 
constraints  which  cannot  be  related  directly  to  the  exit  surface.  In  the 
second  approach,  the  thrust  is  written  in  terms  of  forces  acting  along 
the  wall  contour  AC  in  Figure  1  and  the  entire  region  ABC  is  considered 
in  the  problem.  This  approach  permits  the  consideration  of  flows  with 
dissipative  effects  such  as  gas-particle  drag  or  finite  rate  chemistry. 
This  second  approach  will  be  employed  in  this  work. 

2.  THE  ROW  MODEL 

The  governing  equations  for  the  axisymmetric  gas-particle  flow 
analysis  are  given  in  Ref. (15).  These  equations  are  a  gas  continuity 
equation,  two  system  momentum  equations,  a  system  energy  equation,  a 
particle  continuity  equation,  two  particle  momentum  equations  and  a 
particle  energy  equation. 


(ypu)x  +  (ypv)y  =  0 


0) 


p(uux  +  vuy)  +  px  +  APp(u  -  Up)  =  0 


(2) 
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(3) 


0{uvy  +  VVy)  +  py  +  APp(V  -  Vp)  =  0 

upx  +  vpy  "  a^upx  +  Vpy)  “  =  0  M 

HVx  +  (y»pvp)y  =  o  (5) 

epCup<up>x  +  vp(up>y  -  A<“  '  up»  *  0  <5> 

optyy* +  yy,  -  a<v  -  vp)3  =  °  <?> 

yy  yx +  yyy  - §■  ac<t  -  y  ]  *  °  <«> 

The  parameters  A,  B  and  C  were  defined  for  convenience  and  represent 
particle  drag  ana  energy  parameters. 

A  =  (9fy)/(2mpr2)  (9) 

B  =  (y  -  l)[(u  -  up)2  +  (v  -  vp)2  -  |  C(T  -  Tp)]  (10) 

C  =  (gcp)/(f  Pr)  (11) 

where 

f  *  c(/^cD^Stokes  *12’ 

3  -  N“/("u>stokeS  (,3) 


Assumptions  are  made  that  the  parameters  A  and  C  are  constant  at  least 

2 

locally  in  the  flow  field,  and  relationships  defining  T,  T  and  a  are 

P 

T  =  p/pR  (14) 

tp  ■  +  rc  %  -  0  <15> 

a2  =  2H- 

P 


06) 


hp  and  represent  the  reference  enthalpy  and  temperature  respectively, 
and  Cc  is  the  specific  heat  of  the  condensed  phase.  Cc  depends  on  the 
phase  and  is  equal  to  for  liquid  particles,  «  during  phase  changes, 
and  Cy  for  solid  particles. 

Using  equations  (9)  through  (16),  the  eignt  governing  equations, 

(1 )  through  (8),  can  be  expressed  in  terms  of  two  independent  and  eight' 
dependent  variables.  This  suggests  the  use  of  the  Method  of  Character¬ 
istic!,  for  the  solution  of  the  flow  field  problem.  The  following 
characteristic  and  compatibility  equations  are  obtained  when  the  Method 
of  Characteristics  is  applied  to  the  eight  governing  equations  (1) 
through  (8). 


Along  the  gas  streamline 


dy  _  v 
dx  ~  u 


pudu  'r  pvdv  +  dp  =  -  Ap  [(u  -  u  )dx  +  (v  -  v  )dy] 

P  P  P 


udp  -  a  udp  =  ABppdx 


Along  Mach  lines 


(17) 

08) 

(19) 


3^  =  tan(e  +  a) 


2,  c 

a  (vdu  -  udv)  ±  —  cotodp  -  (udy  -  vdx)  --  - 


-  a-£ 


Along  Particle  Streamlines 


dx  u 

P 


(20) 


|B(udy  -  vdx)  +  a2[(u  -  u  )dy  -  (v  -  v  )dx]l  (21) 

V  r  »  r  J 


(22) 
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(23) 


u„du„  =  A(u  -  u_)dx 

PH  ? 


!pdvp  =  A(v  -  Vp)dx 

(24) 

,pdhp  =  |  AC(T  -  Tp)dx 

(25) 

where  9  is  the  flow  angle  and  a  is  tne  Mach  angle.  The  upper  signs  in 
equations  (20)  and  (21)  refer  to  right-running  Mach  lines  and  the  lower 
signs  refer  to  left-running  Mach  lines. 

This  system  of  equations  provides  seven  equations  along  four  dis¬ 
tinct  characteristic  directions,  and  thus  cannot  be  used  alone  to  solve 
for  eight  variables.  The  reason  for  the  deficiency  has  been  explained 
by  Sauerwein  and  Fendell  (16)  as  being  the  assumption  that  particles  do 
not  contribute  to  the  pressure.  This  results  in  only  one  distinct 
characteristic  direction  for  the  particle  equations,  the  direction  of 
the  particle  streamline.  However,  it  can  be  seen  that  the  particle 
density  term  in  the  particle  continuity  equation  is  also  dependent  on 
the  divergence  of  s  eamlines.  The  absense  of  a  pressure  term  prevents 
expression  of  the  divergence  in  a  characteristic  set  of  equations. 

Further  examination  of  the  system  of  equations  (17)  through  (25) 
reveals  that  if  the  particle  density  can  be  determined  by  other  means,  the 
number  of  dependent  variables  is  reduced  to  seven,  and  the  system  of 
equations  is  adequate.  For  this  purpose,  a  stream  function  is  introduced. 
Since  the  total  amount  of  particles  passing  between  the  centerline  and  a 
given  particle  streamline  must  remain  constant  and  the  distribution  of 
particles  passing  the  initial-value  line  can  be  determined,  an  integra¬ 
tion  technique  can  be  applied  to  determine  particle  density  at  each 


* 
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point  in  the  flow  field.  With  the  particle  density  determined,  the 
remaining  seven  variables  can  be  determined  by  application  of  the 
characteristic  and  compatibility  equations  (17)  through  (25). 

In  the  above  technique,  it  is  necessary  to  determine  the  particle 
streamline  to  find  the  particle  density.  However,  the  particle  stream¬ 
line  depends  on  particle  velocities,  which  have  not  yet  been  determined. 
Thus,  it  is  necessary  to  iterate  the  solution  at  each  point  until  a 
sufficient  accuracy  is  achieved. 

3.  FORMULATION  OF  THE  OPTIMIZATION  PROBLEM 

In  the  problem  of  obtaining  a  maximum  thrust  nozzle,  only  changes  in 
the  nozzle  wall  between  points  A  and  C  (Fig.  2)  in  the  supersonic  region 
are  considered.  The  contour  of  the  nozzle  inlet  and  throat  region  are 
assumed  to  be  fixed.  The  boundary  of  the  problem  consists  of  three 
line  segments.  The  nozzle  wall  defines  the  first,  AC.  Since  the  nozzle 
contour  upstream  of  A  is  fixed,  the  flow  field  upstream  of  a  right-running 
characteristic  attached  at  A  is  independent  of  the  optimization  problem, 
and  the  boundary  AB  is  defined  along  a  right-running  characteristic. 

The  third  line  segment,  BC,  closer  tne  region  and  intersects  the  line  AC 
at  the  end  of  the  nozzle.  No  other  restrictions  are  made  on  this  line  in 
the  formulation  of  the  problem.  However,  the  line  BC  is  represented  as 
a  left-running  characteristic  In  Figure  2  because  the  solution  of  the  cal¬ 
culus  of  variations  equations  results  in  the  characteristic  equation  for 
this  line,  as  shown  later  in  equation  (64). 

The  particles  flowing  through  the  nozzle  do  not  follow  the  wall 
contour,  and  there  will  be  a  region  near  the  wall  which  contains  only 
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gas.  The  limiting  particle  streamline  is  represented  in  Figure  2  by  the 
line  DE.  By  definition  the  nozzle  is  not  optimum  if  the  limiting  parti¬ 
cle  streamline  impinges  on  the  nozzle  wall  since  the  nozzle  would  be 
physically  damaged,  and  in  the  case  where  the  maximum  thrust  solution 
tends  to  result  in  impingement,  it  is  necessary  to  change  the  inlet  and 
throat  region  geometry  sufficiently  to  provide  the  necessary  separation. 

The  thrust  developed  along  the  nozzle  wall  is  given  by 


where  n  is  the  radial  coordinate  of  the  wall  and  is  a  function  of  x. 

This  quantity  is  10  be  maximized.  However,  there  are  several  restric¬ 
tions  which  must  be  imposed  on  the  solution.  The  first  requires  that  the 
wall  be  a  gas  streamline.  Thus, 

un'  -  v  =  0  on  AC  (27) 

Equation  (27)  is  multiplied  by  np  for  convenience  in  subsequent  analysis 
to  give  the  equivalent  expression 

np(un‘  -  v)  =  0  on  AC  (28) 

The  second  restriction  is  a  constraint  on  the  nozzle  wall  such  as  a 
constant  length,  a  constant  arc  length,  etc.  For  problems  of  general 
interest,  this  isoperimetric  constraint  can  be  expressed  in  terms  of  n, 
n*  and  p 

C 

j  G(n,n'.p)dx  -  Constant  =  0  on  AC  (29) 

A 

n 


Finally,  the  governing  equations  (1)  through  (8)  are  restrictions  through¬ 
out  the  region.  Using  L.  to  represent,  the  eight  governing  equations,  the 
functional  which  is  to  be  maximized  is 

fC 

1  =  C(p  -  P0)nn'  +  C-jG  +  C2np(un‘  -  v)]dx 
A 

ff  8 

+  I  h-L.  dy  dx  (30) 

"ABC  i  =  l  1  1 

where  C-j  is  a  constant  Lagrange  multiplier,  C2  is  a  Lagrange  multiplier 
which  is  a  function  of  x,  and  h..  ,  (i  =  1,8),  are  Lagrange  multipliers 
which  are  functions  of  x  and  y.  The  constant  in  the  isoperimetric 
constraint  does  not  affect  the  optimization  problem,  and  therefore  has 
been  dropped  from  I.  For  convenience  in  the  algebraic  development,  the 
following  definitions  are  made: 

Hj  =  (p  -  P0)nn'  +  C-jG  +  C2np(un'  *  v)  (31) 

H2  n,Lt  (32) 

4.  CALCULUS  OF  VARIATIONS 

a.  Necessary  Conditions.  Calculus  of  variations  provides  a  set 
of  conditions  which  must  be  satisfied  if  the  optimum  solution  has  been 
achieved.  Necessary  conditions  include  the  Euler  equations  applicable 
to  the  region  of  interest,  the  transversal ity  condition  applicable 
along  the  boundaries,  the  corner  condition  applicable  to  corners  formed 
by  the  intersection  of  two  boundary  line  segments,  and  the  Erdmann-Heier- 
strass  condition  applicable  to  corner  lines  in  region  ABC.  In  a  super¬ 
sonic  nozzle,  physical  considerations  rule  out  the  class  of  solutions 


uhifh  inrliiHA  i>/»>nap  linac*  ^ne  fko  Cp/lma Rn^Uaiopetpa tr  rnnHitian  it 
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not  used  in  this  formulation  and  the  first  three  conditions  are  consid¬ 
ered  sufficient. 

b.  The  Euler  Equations.  The  Euler  equations,  as  given  by  Hiele 
(17),  are 


(k  =  1,2*..., 8) 


(33) 


where  z ^  represents  the  eight  dependent  variables  p,  p,  u,  v,  hp,  p^,  u^, 
and  Vp,  and  pj,  and  q^  represent  the  derivatives  of  these  variables  with 
respect  to  x  and  y  respectively.  The  Euler  equations  are  developed  in 
Appendix  I.  The  results  are 

-  Vx  -  Vx  -  h4  l  <px  -  a2px>  +  *(hl>x  +  tt<h2>x  +  ,(h2>y  = 

*  h2  1+  A  Pf  [h2  -  h42{t-l  )(°-up)  -  Kg]  <*» 

-  -  Vy  '  h4  p  tpy  -  aV  +  y<h,)y  +  “(h3ix  +  v<fc3)y  = 

*  hj  i  +  A  °f  [h3  -  h42(y-l  )(.-«p)  -  h7]  (35) 


h4“x  +  h4,y  +  h4  r  (o<,x  *  v°y)  *  (hJx  +  *h3!y  *  u{fe4)x  + 


♦v'h4)y*A  J.||[(y-l)h4-h8] 


(36) 


rpx  +r‘py  +  3™%>x  -  yv(hl  >y  +  a2<t,2ix  +  aZ(h3>y  * 

=  A  (-rl)CT[(y-l)h4  -  h8]  +  yi4B  -  h2{u-up)  -  h3(v-/p)i  (37) 

h6(up>x  +  h7(vP)x  +  h8(hp!x  *  >(h5>x  -  Wx  -  Vh6>y  = 

=  A(hj  -  h42(y-l  )(u  -  up)  -  h6)  -  h6  ^  (38) 
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h6(up>y  +  VVy  +  Wy  '  ^Jy  '  up(h7>x  '  Vp(h7)y 
.A^3-h42(Y-l)(*-vp)-h7}-h7^ 

Vs’x  +  V^y  ■  !  ^  {"8  •  W  )h4}  +  h8  ^5- 


i“p(h6>x  +  vp<h5>yJ 


■  A{h2(u  '  V  +  h3(v-vp>  "  h48} 


(39) 

(40) 

(41) 


In  equations  (38)  through  (41)  the  particle  density  function,  p  ,  has 

r 

been  factored  out.  It  must  be  remembered  that  these  equations  apply 
only  in  the  region  where  particles  are  present,  region  BDE.  Since  the 

dependent  variables  are  determined  independently  of  the  optimization  prob¬ 
lem,  they  can  be  considered  as  known  values  durina  the  evaluation  of  the 

Euler  equations.  Thus  equations  (34)  through  (41)  constitute  eight 
equations  with  eight  unknown  variables,  h^  through  hg.  Application  of 
the  Method  of  Characteristics  (See  Appendix II)  results  In  a  set  of  seven 
compatibility  equations  applicable  along  four  distinct  characteristic 
directions. 

Along  the  gas  streamline 


&=  1 
dx  u 


(4?.) 


1  2 

-  I^du  -  hgdv  -  h^dp  -  a  dp)  +  ydh-j  +  udhg  + 

+  vdJij  =  |h2  y  +  A^£hj  -  h42(y-l)(u-up)  -  h6]jdx  + 

+  {l>3  j  *  A^{h3  -  h42(yl  )(v-vp)  -  h^ljdy  (43) 


* 


14 


.  ■  f -rJ  'MSSHSHypHB  m  ".II- ,. 1  JRL'M' .'.  -l » u  ■„  .WPajW|PMJ.PS^^WPWll 


-  n3ujav  *  n2  -  ap  -  iY-Uua~n4  -  dp  +  rfuan-|  - 
-  ua2dh4  =  |-a2h4  p-  -  *2A  f- §  |  C(y-1  )h4  -  hQ3  + 

+  A  I  ^Y_1  ^ct^y-1  >h4  -  hg]  +  A  [Yh4B  -  h2(u-up)  -  h3(v-vp)] 

• ft  ^  (h2  a^-- h3)(v-vp)}dx  (44) 

Along  gas  Mach  lines 

g£=tan(e  +  «)  (45) 

h2du  +  h3*iv  +  h4  ^  (dp  -  a2dp)  -  ydh-j  +  tana(vdh2  -  udh3)  = 

=  ±  tana|h3  J  dx  +  A  ^  [h3  -  h42(Y-l )(v-vp)  -  h?]dx  - 

*  h2  y  dy  *  A  ^h2  ‘  ^(T-Wu-Up)  -  hg]dy  + 

+  ^  (udy  -  vdx)A  ^  (y-1  )[|  Ct|(y-1  )h4  -  hfl|  +  h4B]|  (46) 

where  the  upper  signs  in  equations  (45)  and  (46)  are  applicable  along 
right-runnir.g  characteristics  and  the  lower  signs  are  applicable  along 
left-running  character* sties. 

Along  particle  streamlines 

&  =  !& 
dx  up 

yyfo5  =  A[h2(u-up)  +  h3(v-vp)  -  h4B]dx 

updh6  *  Vh7  =  h6dup  +  h7dvp  +  h8dhp  “  ydh5  ’ 

-  A|h2  -  h42(Y-l )(u-up)  -  hgjdx  +  h&  ^  dx  - 

-  A«|h3  -  h42(Y-l)(v-yp)  -  h?|dy  +  h?  ^  dy 


(47) 

(48) 

(49) 
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Updh8  =  Ir-  ^-h8  '  ("HJh^dx  +  h8  /  dx  (5°) 

Equations  (42)  through  (50)  are  a  system  of  seven  compatibility 
equations  which  apply  along  four  distinct  characteristic  directions. 

Since  there  are  eight  unknown  dependent  variables,  h-j  through  hg,  a 
deficiency  exists.  Analysis  of  the  equations  shows  that  the  system  of 
equations  does  not  permit  evaluation  of  the  variables  hg  and  h7  since 
derivatives  of  both  appear  in  equation  (49)  only. 

In  this  analysis,  hg  and  h-  will  both  be  determined  by  finite 
difference  evaluation  in  the  vicinity  of  the  point  being  considered. 

c.  The  Corner  Condition.  The  corner  condition  applicable  at  the 
corners  A,  B  and  C,  as  given  by  Miele  (17),  is 

3H,  3H, 

A[H,  -  y‘  gyrlta  +  A[^j-36y  =  o  (51) 

where  A  denotes  the  difference  between  the  value  of  the  quantity  in 
brackets  evaluated  on  each  side  of  the  corner  and  5x  and  Sy  signify 
small  variations  in  x  and  y. 

At  point  A,  both  x  and  y  are  fixed,  and  6x  and  6y  are  both  zero; 
thus  the  corner  condition  is  satisfied  without  new  conditions.  At  point 
B,  the  quantity  H-j  is  zero  on  both  sides  of  the  corner,  and  no  new 
conditions  result  At  point  C  the  variations  Sx  and  Sy  are  arbitrary. 
Each  must  be  able  to  vary  independent  of  the  other  so  the  coefficients 
must  be  zero.  Since  is  zero  along  BC,  the  coefficients  of  5x  and  iy 
are  zero  on  that  side  of  the  comer,  and  thus  must  also  be  zero  at  C 
when  evaluated  along  AC. 

3H 

*,  -  y  IF  *  0  *‘C  (52) 
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3H 

^r  =  0  at  C  (53) 

Equations  (52)  and  (53),  when  applied  at  point  C  (see  Appendix  III)* 
result  In 


ci  - 


(pc  -  po)nc\ 


u  fi 
c  c 


<pc  -  po>n  •  & 


r  36 


0  u 
Hc  c 


(54) 

(55) 


d.  Thi.-  Transversal ity  Condition.  Along  the  boundaries,  the  trans- 
versality  condition  must  be  satisfied.  The  transversal ity  equations  are 
given  in  Appendix  IV,  and  when  applied  along  BC  result  in  the  set  of 
equations: 


h3  -  y'h2  +  h4(v  -  y'u)  =  0 

(56) 

^ytv  -  y'u)  -  h4a2(v  -  y'u)  =  0 

(57) 

h2(v  -  y'u)  -  ^yy'  =  0 

(58) 

h3(v  -  y'u)  +  h^y  =  0 

(59) 

h5(y'up  -  vp>  =  0 

(60) 

h6(y‘up  -  V  +  h5y'y  =  0 

(61) 

h7(y'Up  -  vp)  -  h5y  =  0 

(62) 

h8(y‘up  '  V  =  ° 

(63) 
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Substituting  aquations  (57)  through  (59)  into  (56),  we  obtain,  after 
simplifi  cation 


«•  - 


vu  ±  a2\4^_-  1 
u"  -  a2 


on  BC 


(64) 


Thus  the  transversal ity  conditions  specify  that  the  line  BC  must  be 
a  left-running  characteristic. 

Since  BC  is  a  characteristic  neither  of  the  quantities  v  -  y'u  or 
y’Up  -  Vp  are  zero.  Equations  (60)  through  (63)  then  reduce  to 


h5  =  0 

(65) 

h6  =  0 

(66) 

h7  =  0 

(67) 

h8  =  o 

(68) 

and  equations  (56)  through  (59),  in  addition  to  specifying  the  characte 
istic  direction,  become 


h^y  +  h2(y'u  -  v)  =  o 

(69) 

h^yy'  -  h3(y'u  -  v)  =  0 

(70) 

V  -  *\  = 0 

(71) 

When  the  ;rans versa 1 ity  condition  is  applied  along  the  contour  CA 
the  following  eouatlons  result: 

h1  =  C2  (72) 
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J7JLJ!!. JW»V-<  J*-*--- '» J  yJ- 


vn0  +  vn  +  uC,G„ 

2  1  p 


0 


\/3) 


^1  _  du  ^1  r/»  d_  /9G  \  _  /9G\ /dv\-| 

3x  3x  npu  ~  dx  'irT*  “  pJ'3p^dx'~ 


(74) 


5.  FINITE  DIFFERENCE  EVALUATION  OF  hg  AND  h? 

The  failure  of  the  method  of  characteristics  to  produce  eight  com¬ 
patibility  equations  necessitates  the  use  of  another  method  to  solve  for 
at  least  one  unknown  dependant  variable.  Since  equation  (49)  includes 
both  hg  and  hy  in  differential  form,  at  least  one  of  these  must  be  obtained 
by  another  method  in  order  to  permit  the  solution  of  the  other.  In  this 
analysis,  both  will  be  obtained  in  the  vicinity  of  a  point  by  using  the 
appropriate  Euler  equations  and  the  definition  of  a  derivative. 

The  two  Euler  equations  which  contain  partial  derivatives  of  hg  and 
h7  are  equations  (38)  and  (39),  repeated  here  for  convenience 


Wx  +  h7(vp>x  +  Wx  -  y(h5)x  -  Vh6>x  -  vp<h6)y  = 

-  A(h2  -  h42(y-l)(u-up)  -  h6)  -  h6  i  (38) 

"6(up>y  +  h7(Vy  +  h8!hp!y  *  *<h5>y  *  "p^x  '  vp(h7,y  = 

•  a(|>3  -  h42(y-l)(v  -  vp)  -  h?}  -  h?  ^  (39) 

These  equations  can  be  written  in  a  total  differential  form  along  the 
particle  streamline  if  all  the  partial  derivatives  except  the  partial  deri¬ 
vatives  of  hc  and  h7  are  known.  To  solve  for  the  differential  (u  1  ,  equa- 

o  /  p  X 

tion  (6)  and  the  total  derivative  of  uQ  are  used. 
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to  obtain 


VVx  *  Wy  =  A(u  '  V 

duP  °  (up)xdx  +  'up)ydy 


A{u  -  u„)dy  -  v  du 

(y  )  =  - E__ _ 2 _ fl 

P  x  u  dy  -  v  dx 


Along  the  particle  streamline,  (u  )  is  indeterminate.  However,  under 

P  * 

the  assumption  of  continuous  derivatives,  (u  }  can  be  evaluated  along  a 
gas  right-  or  left-running  characteristic  and  the  results  employed  in 
equation  (38). 

In  similar  fashion 

A(v  -  V  )dy  -  v  dv„ 

(y  )  =  — . -  P. . -  P _ 2.  t  -j-»' 

v  p'x  u„dy  -  v„dx 
P  P 

4  AC(T  -  T  )dy  -  vdhn 

(hpJx  -  Updy  -  5pdx  P  P  <78> 

(, )  .  -  A-(u.:.>)dx '  “A.  m 

p  y  Up^y  -  vpdx  l#*' 

A(v  -  v„)cix  -  u  dv„ 

*  Vy  =  Updy  -  vpdxP  * 80) 


<Vy  *  - 


f  AC(T  -  T  )dx  -  u  dh 


u  dy  -  v  dx 
P  p 


Using  equation  (41)  and  the  definition  of  the  total  derivative  of  hg, 
the  following  expressions  are  derived: 


<Vx  • 


A[h2(u-Up)  +  h3(v-vp)  -  h4B]dy  -  yvpdh5 


yupdy  -  yVpdx 
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a.ji 


(^r)w 


A[h2(u-up)  +  hg(v-vp)  -  h4i 


Bjdx 


yupdh5 


tiii  /Tt 

J~P~J 


tin  /TtT 

jrTp- 


(83) 


In  equations  {76)  through  (81),  the  evaluation  of  the  partial  deri¬ 
vatives  involves  only  gas  and  particle  properties.  These  values  are 
known  at  all  points  in  the  flow  field  during  evaluation  of  Lagrange 
multipliers,  and  are  constant.  However,  evaluation  of  equations  (82) 
and  (83)  requires  the  use  of  multiplier  values  not  yet  determined. 

Thus,  values  must  ^irst  be  assumed  for  the  multipliers  and  the  solution 
process  iterated  with  updated  values  for  the  multipliers  until  a  suffi¬ 
ciently  accurate  solution  is  obtained. 

Treating  the  partial  derivatives  evaluated  by  equations  (76)  through 
(83)  as  constants,  and  restricting  equations  (38)  and  (39)  to  particle 
streamlines,  we  obtain  the  two  compatibility-like  equations; 


uP  dh6 =  £  -VVx  -  MVx  "  h8( hp>x  +  y(h5>x + 


+  A  {h2  -  h42(y-l)(u-u  )  -  h6>  -  hg-£]dx 


(84) 


up  dh?  =  [-h6(up)y  -  h7(vp)y  -  hg(hp)y  +  y(h5)y  + 

+  A  (h3  -  h42(y-l)(v-vp)  -  h7)  -  h7-^]dx  (85) 

6.  SUMMARY  OF  RESULTING  EQUATIONS 

The  set  of  equations  which  must  be  satisfied  as  necessary  conditions 
for  the  calculus  of  variations  includes  eight  Euler  equations  or  their 
equivalent  characteristic  and  compatibility  system,  three  transversal ity 
equations  along  the  nozzle  contour  AC,  seven  transversal ity  equations 
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along  the  exit  characteristic  BC,  and  two  corner  condition  equations  at 
the  point  C.  This  is  an  overdetermi ned  system  of  equations.  One  of  the 
transversal ity  equations  along  AC  is  nov  used  in  the  computation  of  the 
Lagrange  multipliers,  and  will  thus  be  available  for  determination  of  an 
error  function.  This  error  function  is  used  to  modify  the  nozzle  contour. 

The  resulting  equations,  in  approximate  order  of  usage  for  computa¬ 
tion  of  multipliers,  are  repeated  here.  At  the  corner  C 


v  % 

c  c 


iPC-p„>D  -vyf.>c] 


Along  the  exit  characteristic  BC 


hly  +  My'u  “  =  0 


hlyy‘  "  h3(y‘u  -  v)  =  0 


V  -  a  =  0 


h5  =  h6  •  h,  =  h8  =  ° 


(65-68) 


Along  the  nozzle  wall  AC 


*’l  _  du  (pc  po)ncvc  1  rr  d  /36  \  3G-* dv-4  i-i»\ 

GT  H5T  - - uX - pun  [6n  ‘  ~  pu(3p'^j  (74) 


Along  gas  Mach  lines 


=  tan (9  +  a) 


* 


h2du  +  h3dv  +  h4  ^  (dP  “  a*dp)  *  yd*^  +  tana(vdh2  -  udh3)  = 

=  +  tana|h3  y  dx  +  A  ^  [h3  -  h42(y-l)(v-Vp)  -  h?]dx  - 

-  h2  ydy  -  A  ^  [h2  -  h42(Y-D(u-up)  -  hgjdy  + 

+  !*■  (udy  -  vdx)A  (y-l)[|  CT  j(y-l)h4  -  hg|  +  h4B]|  (46) 

a  p  '  } 

Along  the  gas  streamline 

&  =  v  (42) 

dx  u 

-  h2du  -  hgdv  -  h4  ^  (dp  -  a2dp)  +  ydh-j  +  udh2  + 

+  vdh3  =  |h2  —  +  A  -jj  [h2  -  h42(y-l )(u-Up)  -hg]j  dx  + 

+  |h3  y  +  A  ^  [h3  -  h42  (y-l)(v-vp)  -  h-]j  dy  (43) 

1  2  1 

(h2v  -  h3u  )dv  +  h2  -  dp  -  (y-l)ua  h4  g-  dp  +  yudh^  - 

-  ua2dh4  =  |  -  a2h4  y  -  a2A  |  §  [(y-l)h4  -  hg]  + 

+  A  ^  |  (-H)  CT[(r-!)h4  -  h8]  +  A  ^  hh4B  -  h2(u-up)  - 

-  h3(v  -  vp)]  -  A  ^  (h2  ^  -  h3)(»  -  vp)|  dx  (44) 

Along  particle  streamlines 

&  =  ^  (47) 

W  up 

yUpdhg  =  A[h2(u-Up)  +  h3{v-vp)  -  h43]dx  (48) 
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(84) 


updhS  =  t'h6(Vx  ‘  h7(vp>x  •  h8(hp>x  +  *(tl5!x  + 

+  A  (h2  -  h42(-H)(u-up)  -  h6)  '  h6  y^* 

Updh?  =  [-h6(up)y  -  h7(vp)y  -  h8(hp)y  +  y(h5)y  + 

+  A  (h3  -  h42(Y-l)(v-vp)  -  h?}  -  h?  ^]dx  (85) 

updh8aft;  ["S'  b-Dh4]dx  +  h8^ix  (50) 

Finally,  equation  (73),  which  is  not  employed  in  the  above  set,  is  avail¬ 
able  as  a  check  condition  applicable  along  the  nozzle  wall  to  determine 
if  the  nozzle  is  optimum. 


vh0  +  vn  +  uC,G„  =  0 
Z  1  p 


(73) 
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SECTION  III 
NUMERICAL  METHODS 


1.  SOLUTION  PROCEDURE 

The  solution  procedure  consists  of  estimating  a  nozzle  contour,  per¬ 
forming  a  flow  field  analysis,  computing  Lagrange  multipliers,  evaluating 
the  check  condition  given  by  equation  (73),  and  modifying  the  nozzle. 

This  procedure  is  repeated  using  the  modified  nozzle  contour  as  the 
estimated  nozzle  contour  until  the  error  criterion  is  satisfied.  This 
section  describes  the  program  developed  for  use  in  the  parametric  studies. 

The  initial  nozzle  contour  may  be  input  in  tabular  form  or  calculated 
internally  in  the  form  of  conical,  parabolic  or  circular  arc  nozzles.  The 
initial-value  line,  required  to  begin  the  supersonic  flow  field  calcula¬ 
tion,  may  be  input  as  data  or,  for  standard  ccnverging-di verging  nozzles, 
can  be  generated  internally  based  on  tha  combustion  chamber  conditions. 

2.  aOW  FIELD  ANALYSIS 

The  flow  field  analysis  utilizes  the  method  of  characteristics. 

Mesh  points  are  located  at  intersections  of  left-running  characteristics 
(LRCs)  and  right-running  characteristics  (RRCs).  The  reference  stream¬ 
line  points  are  obtained  by  linear  interpolation  between  reference  points 
along  the  LRC  and  the  RRC.  See  Figure  3.  The  mesh  construction  proceeds 
along  LRCs  beginning  at  either  an  initial-value  point  or  an  axis  point, 
inserting  points  along  RRCs  from  points  cn  the  previous  LRC  until  the  LRC 
intersects  with  the  nozzle  contour.  Additional  LRCs  are  constructed  in 
this  manner  until  the  end  of  the  nozzle  contour  is  reached. 
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Axis  points  are  located  at  the  intersection  of  an  RRC  from  the  second 
point  on  the  previous  LRC  and  the  axis.  A  limiting  particle  streamline 
point  is  calculated  along  each  LRC  and  at  least  one  additional  point  must 
lie  between  the  limiting  particle  streamline  and  the  nozzle  contour. 

Point  insertion  routines  are  used  to  insert  additional  mesh  points  when¬ 
ever  the  mesh  size  exceeds  the  values  specified  by  input  parameters. 
Whenever  «  mesh  point  is  inserted  along  an  RRC,  a  new  LRC  is  begun  at  that 
point  and  proceeds  to  the  nozzle  contour  before  continuing  on  with  the 
interrupted  LRC.  The  point  insertion  routine  is  also  used  to  locate  the 
last  point  at  the  end  of  the  nozzle  contour. 

During  the  flow  field  analysis,  a  secondary  start  line  is  constructed 
which  follows  an  RRC  from  a  point  on  the  nozzle  contour  near  point  A 
(Fig.  2)  to  the  axis.  This  secondary  start  line  is  used  in  subsequent 
iterations  to  reduce  the  calculation  time  by  eliminating  the  recalculation 
of  mesh  points  unaffected  by  the  wall  modification. 

3.  OPTIMIZATION  CALCULATIONS 

The  optimization  portion  of  the  program  begins  with  the  determination 
of  properties  along  an  LRC  which  passes  through  the  end  point  of  the 
nozzle  contour.  This  is  accomplished  by  linear  interpolation  along  the 
RRCs  which  connect  points  on  either  side  of  the  exit  characteristic.  The 
Lagrange  multipliers  are  then  calculated  along  the  exit  characteristic 
and  the  method  of  characteristics  is  employed  to  determine  the  Lagrange 
multipliers  throughout  the  flow  field  after  reassembling  the  mesh  con¬ 
struction  in  reverse  order.  A  check  condition,  determined  from  equa¬ 
tion  (73)  and  expressed  in  terms  of  slopes,  is  calculated  at  wall 
points  and  an  error  function  is  determined. 


4.  RELAXATION  TECHNIQUE 

The  check  condition,  used  to  determine  if  the  nozzle  contour  is  op¬ 
timum,  will  not  in  general  be  satisfied.  After  rewriting  equation  (73), 


an  error  parameter  is  defined: 


E  =  tan'Un')  -  tan 


-1  h3+Cl6! 


(h2  *  n  J 


By  assuming  that  the  remaining  variables  do  not  change  appreciably  with 
changes  in  wall  slope,ri',  the  error  parameter  may  be  relaxed  by  calculat¬ 
ing  a  new  slope  so  that  E  will  be  zero,  thus,  the  new  slope  is  given  by 

h,+C,Gn 

-  .A  L  P.  /ft R) 


This  relaxation  technique,  first  proposed  by  Major  A.  A.  Taylor,  USAF, 

while  he  was  a  graduate  student  at  Purdue  University,  permits  a  simple 

integration  to  determine  the  relaxed  nozzle  contour. 

* 

n(x)  =  n(xft)  +  n’(x)dx  (89 

■*XA 

Since  the  other  variables  do  change  when  n'  is  changed,  the  solution  is 
not  exact,  and  the  flow  and  optimization  calculations  must  be  iterated. 
It  was  found  that  the  relaxation  scheme  resulted  in  an  overcorrec¬ 
tion  in  the  majority  of  cases.  For  this  reason,  a  weighing  factor  was 
introduced  to  reduce  the  number  of  iterations,  and  thus  the  computation 
time.  Thus,  the  modified  contour  slope  n‘  is  given  by 


q'  =  rt]  +  0.8(rj^  -  o-j)  (9 

where  n),  represents  the  slope  found  by  equation  (88)  and  nj  represents 
the  slope  of  the  original  estimated  nozzle  contour.  In  addition,  res¬ 
trictions  were  placed  on  the  calculation  of  the  n ew  nozzle  contour  to 
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prevent  a  change  in  slope  greater  than  10  degrees  in  either  direction 
and  to  prevent  negative  slopes  from  being  used. 

The  nozzle  contour  could  become  discontinuous  if  point  A  were  fixed. 
This  is  avoided  by  permitting  the  point  A  to  move  along  the  throat  radius 
of  curvature  until  the  slope  at  point  A  equals  the  slope  calculated  by 
equation  (88).  In  the  event  point  A  moves  upstream  of  the  first  point 
of  the  secondary  start  line,  the  flow  field  analysis  will  begin  at  the 
initial-value  line  for  the  next  iteration  and  a  new  secondary  start  line 
beginning  upstream  of  the  new  point  A  will  be  calculated. 

Convergence  of  the  nozzle  contour  is  determined  by  the  value  of  the 
error  parameter,  E,  determined  by  equation  (87).  When  E  is  less  than  the 
value  of  the  input  parameter  TOL  at  each  point  along  the  nozzle  contour, 
the  convergence  criterion  is  met  and  the  program  is  terminated.  If  this 
convergence  criterion  is  not  met,  the  estimated  nozzle  contour  is  re¬ 
placed  by  the  modified  nozzle  contour  and  the  solution  is  iterated 
beginning  at  the  flow  field  analysis. 

Figures  4  and  5  illustrate  the  behavior  of  this  relaxation  scheme. 
The  sample  case  consists  of  a  design  to  produce  a  maximum  thrust  nozzle 
9.6  inches  long  with  a  throat  radius  of  1.2  inches  and  a  throat  radius 
of  curvature  of  2.4  inches.  The  initial  estimate  of  the  contour  is  a 
25°  conical  nozzle.  Aluminum  oxide  particles  4  microns  in  diameter  with 
a  mass  flow  rate  0.4  times  the  gas  flow  rate  are  specified.  The  gas 
properties  are:  y  -  1.28,  molecular  weight  *  17.76,  chamber  pressure  = 
500  psias  chamber  temperature  =  6500°R,  and  ambient  pressure  =  3  psia. 

The  convergence  tolerance,  TOL,  was  0.05  degrees. 

Seven  iterations  are  required  for  convergence.  The  initial,  second, 
fourth,  and  final  nozzle  contour  are  shown  in  Figure  4.  The  fourth  and 
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final  contours  are  almost  identical.  The  corresponding  error  values  are 
shown  in  Figure  5.  The  initial  contour  estimate  of  a  conical  nozzle 
resulted  in  large  error  values  in  the  first  iteration.  However,  the 
changes  in  nozzle  contour  were  limited  to  ten  degrees.  The  error  values 
in  the  second  iteration  begin  to  show  the  overcorrecticn,  then  reverse 
where  the  first  iteration  correction  was  limited  to  ten  degrees.  The 
error  values  in  the  fourth  iteration  show  a  marked  decrease,  and  by  the 
seventh  iteration  are  less  than  0.05  degrees  at  all  points. 
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SECTION  IV 


PARAMETRIC  STUDIES 

1 .  GENERAL 

A  computer  program,  written  in  FORTRAN  IV  for  the  CDC  6500,  was 
developed  based  on  the  analysis  presented  in  Section  II  and  the  numerical 
methods  presented  in  Section  III.  This  computer  program  was  designed  to 
oermit  either  the  analysis  of  a  given  nozzle  or  the  design  of  a  maximum 
thrust  nozzle  subject  to  some  geometric  design  constraint.  The  program 
user  may  choose  to  begin  at  an  initial -value  line  in  the  supersonic 
region  or  with  combustion  chamber  conditions.  Output  options  permit 
selection  of  the  amount  of  printed  output  and  a  punch  capability  designed 
to  permit  input  of  the  computed  nozzle  contours  and/or  start  lines  into 
subsequent  computations.  A  program  description,  description  of  input 
variables, and  lifting  of  several  sample  rases  are  given  in  Reference  (18). 

This  section  presents  the  results  of  parametric  studies  to  determine 
the  effect  of  varying  several  of  the  input  parameters.  For  each  parametric 
study,  one  variable  was  varied  holding  all  others  constant. 

2.  VARIATION  OF  MESH  SIZE 

The  method  of  characteristics  mesh  construction  is  controlled  in  two 
ways.  First,  the  number  of  initial -line  points  (NILF)  determines  the 
size  of  the  mesh  near  the  initial  value  line  and  affects  the  mesh  size 
throughout  the  flow  field.  Second,  provisions  are  incorporated  into  the 
design  program  to  limit  the  mesh  size.  In  supersonic  divergent  nozzle 
flow,  the  mesh  size  tends  to  increase  as  the  mesh  construction  proceeds 
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by  the  program  whenever  the  dimension  along  a  riqht-running  characteris¬ 
tic  exceeds  the  input  value  DR,  the  dimension  along  a  left-running 
characteristic  exceeds  the  input  value  DL,  or  the  flow  angle  change  along 
the  nozzle  contour  exceeds  the  input  variable  DTWI. 

For  the  first  part  of  this  parametric  study,  the  values  of  DL,  DR, 
and  DTWI  were  set  to  relatively  large  values,  thus  inhibiting  the  mesh 
point  insertion  routines.  A  nozzle  with  a  1.2  inch  throat  radius, 
length  of  9.6  inches  and  a  25  degree  cone  was  selected.  The  ratio  of 
the  mass  flow  rate  of  the  particles  to  the  mass  rate  of  flow  of  the  gas 
was  0.4  ar.d  the  size  of  the  particles  was  4  microns.  Three  values  of 
NILP  were  used  to  vary  the  mesh  size.  The  results  are  shown  in  Table  I. 
The  times  required  for  each  design  were  342.9  sec,  593.8  sec,  and  826.1 
sec  for  NILP  =  12,  NILP  =  16,  and  NILP  =  20  respectively.  These  times 
could  be  reduced  significantly  if  better  initial  contours  were  used.  In 
all  nhree  cases,  the  initial  contour  is  the  same.  The  final  nozzle  con¬ 
tour  was  reached  in  one  less  iteration  for  NILP  =  20. 

Increasing  NILP,  and  thus  decreasing  mesh  size,  results  in  almost 
identical  final  nozzle  contours,  with  the  nozzle  exit  radius  slightly 
larger  for  increased  NILP.  The  values  for  computed  thrusts  indicate  a 
dependence  on  the  mesh  size.  This  can  readily  be  seen  by  comparing  the 
thrust  computed  for  the  first  contour  estimate  in  each  case.  Since  these 
contours  are  identical,  the  difference  can  be  attributed  to  the  effect 
of  the  mesh  size  on  the  numeric  scheme.  Another  effect  can  be  noted  in 
reviewing  the  thrust  values.  The  final  calculated  thrust  is  typically 
slightly  less  than  one  or  more  of  the  intermediate  calculated  thrusts. 
However,  this  effect  is  slight,  amounting  to  approximately  Q.03%  of  the 
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suggests  that  the  tolerance  may  be  increased  so  that  fewer  iterations 
are  required,  since  intermediate  contours  yield  essentially  the  same 
thrust  as  the  final. 

A  second  part  of  this  parametric  study  consisted  of  reducing  the 
calculation  times.  For  this  purpose,  four,  five  ar.d  six  iterations 
(input  parameter  NITER)  were  performed  with  NILP  =  12.  Then  the  result¬ 
ing  computed  nozzle  contour  was  input  into  a  design  program  with 
NILP  =  20,  and  run  until  the  error  tolerance  equal  to  0,05°  was  satisfied. 
The  results,  shown  in  Table II,  indicate  a  savings  in  total  computational 
time  can  be  achieved  by  first  obtaining  a  better  first  estimate  of  the 
nozzle  contour  using  a  large  mesh  size.  However,  increasing  the  nunber 
of  iterations  with  NILP  =  12  to  5  and  6  does  not  reduce  the  nunber  of 
iterations  required  with  NILP  »  20.  A  similar  savings  can  be  achieved 
by  selecting  an  appropriate  parabolic  nozzle  contour  for  the  first 
estimate. 


TABLE  II 

COMPUTATION  TIMES 


NUMBER  OF  ITERATIONS  WITH  NILP  =  12 


NILP 

0 

4 

5 

6 

12 

178.6 

217.3 

253.6 

20 

826.1 

360.6 

361.6 

359.7 

826.1 

539.2 

578.9 

613.3 

3.  VARIATION  OF  PARTICLE  SIZE 


For  the  particle  size  studies,  an  eight  inch  circular  arc  nozzle 
with  a  30°  attachment  angle  and  a  four  inch  exit  radius  was  used.  The 
inlet  to  the  throat  was  30°.  Particle  sizes  of  2,  4,  6,  8  and  10  microns 
were  used.  In  all  cases,  the  ratio  of  the  mass  rate  of  flow  of  the 
particles  to  the  mass  rate  of  flew  of  the  gas  was  0.4.  The  nixnber  of 
initial  line  points  used  was  20. 

The  final  contours  are  shown  in  Figure  6.  The  optimised  contours 
for  the  4,  6,  8  and  10  micron  particles  are  nearly  the  same,  and  the 
2  micron  optimized  contour  is  not  consistent  with  the  others.  Analysis 
of  the  numeric  results  revealed  that  the  limiting  particle  streamline  for 
the  2  micron  case  was  near  the  nozzle  contour,  causing  many  linear  inter¬ 
polations  in  the  process  of  computing  limiting  particle  streamline  points 
and  wall  points.  This  caused  a  loss  of  accuracy. 

An  analysis  of  off  design  conditions  was  made  to  determine  the 
effect  of  particle  size  mismatch.  The  contours  used  in  the  analysis  are 
the  final  contours  just  mentioned.  In  addition  lu  cne  2,  4,  6,  8  and  10 
micron  sizes,  an  analysis  was  also  made  using  equal  parts  by  mass  of  2, 

4,  6  and  8  micron  particles  such  that  the  total  ratio  of  mass  flow  rates 
was  the  same  as  for  the  analysis  with  one  particle  size.  The  results  of 
these  analyses  are  given  in  Table  III.  Also  given  are  the  nozzle  radii 
at  the  exit  for  the  optiiraxn  contours. 

The  use  of  the  2  micron  wall  resulted  in  lower  thrust  than  for  the 
other  walls  regardless  of  the  size  of  particle  used.  This  confirms  the 
conclusion  that  the  2  micron  design  did  not  provide  an  optimum  contour. 
The  4,  6,  8  and  10  micron  walls  provide  the  same  thrust  for  each  particle 
size  used,  within  the  accuracy  of  the  numeric  scheme  of  the  program.  The 
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parametric  study  required  an  inlet  angle  in  the  neighborhood  of  30°  as  a 
compromise  between  the  2  micron  and  10  micron  requirements.  Smaller  in¬ 
let  angles  resulted  in  impingement  of  the  2  micron  particles  on  the  wall 
while  larger  angles  provided  large  vertical  momentum  to  the  10  micron 
particles,  causing  the  particles  to  cross  the  centerline.  The  numeric 
integration  of  the  particle  mass  could  not  be  performed  under  these  cir¬ 
cumstances.  Thus,  in  normal  design  programs  for  small  particle  sizes, 
an  inlet  angle  of  45°  is  recommended  to  provide  a  greater  separation 
between  the  nozzle  contour  and  the  limiting  particle  streamline. 

The  use  of  larger  particle  sizes  results  in  less  thrust,  reflecting 
greater  dissipative  effects.  The  higher  thrust  computed  for  the  mixture 
of  four  particle  sizes  is  due  in  part  to  the  use  of  some  smaller  particle 
sizes,  and  in  part  the  increase  in  mesh  points  caused  by  insertion  of  a 
mesh  point  at  each  limiting  particle  streamline.  The  same  effect  was 
noted  in  the  mesh  size  parametric  studies. 

4.  VARIATION  OF  PARTICLE  MASS  FLOW  RATE 

In  this  parametric  study,  the  ratio  of  the  mass  rate  of  flow  of  the 
particles  to  the  mass  rate  of  flow  of  the  gas  (WPWGT)  was  set  at  0.1, 

0.4,  0.7,  and  1.0.  For  the  larger  value  of  WPWGT,  the  effects  of  the 
particles  near  the  throat  region  result  in  moving  the  sonic  line  down¬ 
stream.  For  this  reason,  the  initial-value  line  was  shifted  downstream 
by  adjusting  the  input  variables  THIW  ana  ZAX.  The  nozzle  is  an  8  inch 
nozzle  with  the  initial  contour  assumed  to  be  a  circular  arc  attached 
at  a  30°  angle  and  with  an  exit  radius  of  4  inches.  The  particle  size 
used  was  4  microns.  NTLP  was  12 

The  results,  shown  in  Figure  7,  indicate  that  greater  amounts  of 
particles  require  more  expansion  to  maximize  tne  thrust.  This  phenomenon 
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VARIATION  OF  PARTICLE  MASS  TO  GAS  MASS  RATIO 


is  uut-  to  a  combination  of  drag  and  heat  transfer  effects:  The  transfer 
of  energy  between  the  particles  and  the  gas  also  results  in  increased 
thrust.  This  study  did  not  consider  changes  in  the  combustion  properties 
due  to  change  in  composition. 

5.  VARIATION  OF  INLET  ANGLE 


When  the  inlet  angle  was  changed,  using  the  inlet  angles  of  12°, 

30°  and  48°,  it  was  found  that  the  optimum  nozzle  contour  was  the  same 
for  all  three  cases.  However,  the  limiting  particle  streamline  was  alter¬ 
ed  as  shown  in  Figure  8,  and  the  thrust  value  was  affected  by  the  oppor¬ 
tunity  for  energy  exchange  between  the  particles  and  the  gas.  A  9.6 
’nch  nozzle  was  specified,  again  using  4  micron  particles  and  WPWGT  =  0.4. 
The  thrust  increased  as  the  particles  occupied  greater  portions  of  the 
nozzles,  permitting  more  effective  energy  transfer.  This  indicates  that 
the  nozzle  inlet  contour  should  be  designed  to  bring  the  limiting  parti¬ 
cle  streamline  near  the  end  of  the  nozzle  contour  while  avoiding  impinge¬ 
ment. 


6.  VARIATION  OF  THE  DRAG  AND  HEAT  TRANSFER  COEFFICIENTS 

Parametric  studies  wer»  conducted  to  determine  the  effect  cf  inac¬ 
curacies  in  the  empirical  drag  and  heat  transfer  coefficients  used  in 
the  program.  The  parametric  studies  wore  performed  using  a  nozzle  length 
of  9.6  inches,  WPWGT  =  0.4,  4  micron  particle  size  and  NILP  =  12.  The 
inlet  angle  was  30  degrees. 

The  drag  coefficients  and  heat  transfer  coefficients  are  calculated 
from  tabular  data.  The  tabular  data  give  empirical  ratios  of  the  coef¬ 
ficients  for  specific  values  of  Reynolds  numbers  to  the  coefficients  for 
Stokes  flow  regime.  The  drag  coefficient  and  heat  transfer  coefficient 
tables  used  in  the  program  are  the  ones  used  by  Kliegel  and  Nickerson  (14), 
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but  can  be  changed  easily  by  the  program  user  if  desired.  Rarefaction 
corrections  are  applied  to  the  drag  and  heat  transfer  coefficients 
calculated  from  the  tables. 

Figure  9  compares  the  effect  of  increasing  the  drag  coefficient  by 
a  factor  of  3  and  decreasing  the  drag  coefficient  by  a  factor  of  3  rela¬ 
tive  to  the  standard  case.  The  optimized  nozzle  contours  were  the  same 
for  all  three  cases,  but  the  limiting  particle  streamline  was  affected 
substantially.  The  case  with  the  higher  drag  coefficient  results  in 
greater  thrust.  Since  the  particles  occupy  a  greater  portion  of  the 
no_zle,  the  heat  transfer  is  more  effective. 

figure  10  compares  the  affect  of  varying  the  heat  transfer  coefficient. 
The  larger  heat  transfer  coefficient  results  in  larger  computed  thrust 
values  and  nozzle  contours  with  greater  expansion  ratios.  A  ninefold  in¬ 
crease  in  heat  transfer  coefficient  increases  the  thrust  0.5  percent. 
Particle  impingement  is  not  of  concern  since  the  limiting  particle 
streamline  moves  in  the  same  direction  and  by  approximately  the  same 
distance  as  the  nozzle  Contour  moves. 

7.  VARIATION  OF  THE  THROAT  RADIUS  OF  CURVATURE 

The  throat  radius  of  curvature  (RRT)  is  expressed  in  throat  radii. 

The  throat  radius  of  curvature  parametric  study  compares  results  with 
throat  radii  of  curvature  of  1.5  throat  radii,  2.5  throat  radii  and  3.5 
throat  radii.  An  inlet  angle  of  30°  and  nozzle  length  of  8  inches  were 
used.  The  computed  thrust  values  are  nearly  the  same  although  the  nozzle 
contours  and  expansion  ratios  vary.  (See  Figure  11),  The  supersonic 
portion  of  the  nozzle  contours  are  very  similar,  however,  if  they  are 
translated  in  the  axial  direction  a  sufficient  distance  to  compensate 
for  the  variation  in  throat  length. 
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FIGURE  9.  E' 


FIGURE  10.  EFFECT  OF  HEAT  TRANSFER  COEFFICIENT  VARIATION 


OF  CURVATURE 


The  limiting  particle  streamline  was  also  affected  by  the  change  in 
the  throat  radius  of  curvature.  With  larger  values  for  RRT,  the  curva¬ 
ture  of  the  wall  at  the  throat  region  is  more  gradual,  permitting  the 
particles  to  accelerate  more  in  the  axial  direction  prior  to  being  sub¬ 
jected  to  vertical  drag  forces.  The  resulting  momentum  provides  a 
greater  separation  between  the  nozzle  wall  and  the  limiting  particle 
streamline. 

This  parametric  study  indicates  that  it  may  be  advantageous  to  use 
the  larger  throat  radius  of  curvature,  since  thrust  is  not  lost,  particle 
impingement  with  the  wall  is  avoided,  and  there  is  a  possibility  that  the 
nozzle  weight  could  be  reduced. 

8.  NOZZLE  SCALING 

Since  two-phase  flow  is  dissipative,  scaling  of  the  results  of  one 
design  to  provide  a  nozzle  of  greater  or  less  thrust  is  not  applicable. 

A  parametric  study  was  conducted  to  illustrate  the  effects  of  scaling. 

An  inlet  angle  of  45°  and  a  throat  radius  of  curvature  of  3  throat  radii 
were  used  for  each  case.  The  nominal  case  was  a  4.8  inch  long  nozzle  with 
throat  radius  of  0.6  inch.  Additional  cases  were  two  times,  four  times 

«Tr 

an*i  eight  times  the  size  of  the  nominal  case  nozzle. 

The  maximum  thrust  nozzle  contours  for  the  larger  nozzles  result 
in  a  greater  expansion  ratio.  (See  Figure  12).  A  more  marked  effect, 
however,  is  seen  in  the  locations  of  the  limiting  particle  streamlines. 

The  location  of  these  streamlines  may  dictate  changes  in  the  fixed  por¬ 
tion  of  the  nozzle  in  order  to  prevent  impingement  of  particles  on  the 
nozzle  contour  or  to  increase  the  region  in  which  energy  can  be  trans¬ 
ferred  from  the  particles  to  the  gas. 
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12.  EFFECT  OF  SCALING 


The  calculated  thrust  for  each  of  the  final  contours  is  shown  in 
Table  IV.  If  scaling  of  the  nozzles  and  the  thrust  were  applicable,  the 
thrust  would  increase  as  the  square  of  the  increase  in  size.  The  re¬ 
sults  show  that  the  thrust  increase  is  slightly  higher  than  direct  scal¬ 
ing  would  indicate. 
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SECTION  V 
CONCLUSIONS 

A  formulation  of  the  maximum  thrust  "isymmetric  gas-particle  nozzle 
problem  was  presented.  A  computer  program  was  developed  based  on  the 
optimization  analysis.  This  computer  program  was  used  to  conduct  several 
parametric  studies  designed  to  determine  the  effect  of  various  design 
parameters  on  the  nozzle  shape,  limiting  particle  streamlines,  and  thrust 
of  optimum  nozzles. 

The  computer  design  program  can  be  used  either  for  analysis  of  nozzle 
flow  or  for  design  of  the  maximum  thrust  nozzle  contour.  Analysis  of 
flow  fields  can  be  performed  with  four  discrete  particle  sizes  while  the 
design  of  maximum  thrust  nozzles  is  limited  to  one  particle  size. 

The  results  of  several  parametric  studies  were  presented.  The  con¬ 
clusions  based  on  these  studies  include: 

1.  The  change  in  nozzle  contours  as  mesh  size  changes  is  negligible 
and  the  change  in  value  of  computed  thrust  is  small.  However,  the  mesh 
size  does  affect  computer  run  time  significantly. 

2.  The  effect  of  particle  size  on  the  final  nozzie  contour  was  in¬ 
significant  except  for  small  particles  whose  limiting  particle  streamline 
lies  near  the  nozzle  contour,  thus  causing  numerical  scheme  errors. 

3.  The  ratio  of  the  mass  rate  of  flow  of  particles  to  the  mass 
rate  of  flow  of  gas  has  a  significant  effect  on  both  the  final  nozzle 
contour  and  the  thrust.  Higher  concentrations  of  particles  require  a 
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greater  expansion  ratio  to  develop  the  maximum  thrust  and  develop  greater 
thrust. 

4.  Variation  of  the  inlet  angle  to  the  nozzle  throat  does  not 
appreciably  affect  the  final  nozzle  contour,  but  doss  affect  the  limiting 
particle  streamline  and  the  thrust.  Smaller  inlet  angles  result  in 
particles  occupying  a  greater  portion  of  the  nozzle  and  a  greater  thrust. 

5.  Increasing  the  values  of  the  drag  coefficient  does  not  appreci¬ 
ably  affect  the  nozzle  contour,  but  increases  the  portion  of  the  nozzle 
which  contains  particles  and  increases  t!,e  rhrust. 

6.  Increasing  the  heat  transfer  coe’.'cierit  increases  the  expan¬ 
sion  ratio  of  the  final  nozzle  contour,  increases  the  portion  of  the 
nozzle  which  contains  particles,  and  increases  thrust. 

7.  Increasing  the  threat  radius  o'f  curvature  decreases  the  nozzle 
diameter  but  dees  not  affect  the  thrust  significantly.  The  increased 
throat  radius  of  curvature  results  in  greater  separation  between  the 
limiting  particle  streamline  and  the  nozzle  contour. 
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APPENDIX  I 

DERIVATION  OF  THE  EULER  EQUATIONS 


The  extremal  problem  for  the  gas-particle  0022!°  optimization  problem 
is  given  by  the  equation 


I  = 


Jvjdx  + 


8 

f 

i=l 


h-L-dydx 


where 


Hj  =  (p  -  PQ)nn'  +  CjS  +  C2np(un'  -  v) 


(  91) 


v  3.) 


\  =  h1[yPux  +  yupx  +  yPvy  +  yvPy  +  Pv]  + 

+  h2Lpuux  +  pvuy  *  px  +  App(u  -  up)3  + 

+  h3[puvx  +  Pvvy  +  py  +  Acp(v  -  vp)]  + 

+  h4[upx  +  vpy  -  a2uox  -  a^vpy  -  ABPp3  + 

+  h5^p(up>x  +  *upVx  *  »p<Vy  +  Wy  +  »pvp]  + 

+  hfiC0pupfup)x  +  Vp(“p>y  -  »pA!“  -  up>j  + 

+  h,[p  u  (v  )  +  p  v  (v  )  -  p  A(v  -  v  )3  + 

7lPp  p'  p'x  Mp  p'  p'y  ^p  v  p/J 

■  h8[Vp(hp,x  +  Vp-Vy  '  I  <%(T  ‘  V5 


The  functional,  I,  reaches  a  maximum  when  the  maximum  thrust  nozzle 
design  has  been  achieved.  Under  these  conditions,  the  functional  remains 
constant  as  minute  changes  are  made  in  each  of  the  dependent  variables 
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possible  then  to  take  variations  by  setting  the  functional  equal  to  the 
same  functional  with  one  of  the  dependent  variables  incremented  an  infin- 
itesmal  amount.  The  Euler  equation  used  in  the  calculus  of  variations 
technique  is  one  method  of  taking  these  variations  ever  the  region  being 
considered  and  results  in  a  set  of  eight  equations. 

The  general  form  of  the  Euler  equation  is  given  by  Miele  (17)  as 


.  3H„  a  3H0 

"  3x  %pj^  '  3y  ^3q^  =  0  =  1,2, ***,8) 


(33) 


where  z.  denotes  the  eight  dependent  variables  u,  v,  p,  p,  u  ,  v.  p^ 

k  p  p  p 

and  hp,  and  pk  and  are  defined  as 


„  _!?» 

pk  3x~ 


(93) 


n  -!fk 

qk  ~  3  y 


(94) 


A  and  C  are  constants,  and 


3  =  (Y-l)[{u-u  )2  ♦  (v  -  v  )2  -  |  C(T  -  T  )] 


T  =  S. 
PR 


T  =  T_°  +  tX  (h  -  h  °) 
P  P  cc  p  p 


(10) 

(14) 

(15) 


,2  =  IE 


(16) 


hp°  and  Tp0  represent  the  reference  enthalpy  and  temperature  respectively 
and  Cc  is  the  specific  heat  of  the  condensed  phase. 


Tne  general  tular  equation,  when  applied  to  equation  i  92)  with 
2^  =  u,  yields 

(hjyp)x  +  h1y&x  +  h2pux  -  {^2oU)x  "  (h2pV^y  + 

+  h2%  +  h3pvx  *  Vx  *  h4a2px  “  Wph-1  Mu  -  up)  - 
'  h6App  =  0  (95) 

Simplifying  and  grouping  terms  results  in 

*  ypfhl)x  +  h2oUx  “  h2^pux  +  upx  +  pvv  +  vpv^  ' 

-  pu(h2)x  -  pv(h2)y  +  h3pvx  +  h4(px  -  a2p  J  + 

+  App[h2  -  h4(y-l)2(u  -  Up)  -  hg]  =  0  (96) 

Using  the  gas  continuity  equation  and  rearranging  gives  the  result 

■  h2Ux  ■  h3Vx  ”  h4  p  {°x  '  a2px)  +  y'Vx  +  u(h2}x  +  v(h2)y  = 

=  h2  y  +  A  )(u*up)  -  h6]  (34) 

Performing  the  same  operations  with  z^  =  v  yields 

-  (hjyp)y  +  h7yPy  +  h,P  +  n2puy  -  (h3pu)JC  + 

+  h3pvy  -  <h3pv)y  +  h3App  +  h4Py  "  V%r  - 

-  h4App(Y-l)(v  -  vp)  -  h7App  -  0  (97) 

-  yp(hj  )y  +  h2p*jy  -  h3[cux  +  upx  +  pvy  +  vPy]  + 

+  h3pvy  -  pu^3}x  ~  py{h3}y  +  h4^y  *  »%>  + 

+  ^PpC^  -  h4(y-l)2(v  -  vp)  -  h?]  =  0  (98) 
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(35) 


h2uy  '  h3Vy  '  h4  p  (py  '  3%}  +  y(hl'y  +  u{h3}x  +  v(h3}y  = 
h3  y  +  A  “p  'h3  -  h42<Y-l)(v-vp}  -  h?] 


For  Zt,  =  p,  the  following  result  is  obtained. 

IV 


(h2Jx  '  <h3>y  -  (h4u)x  -  !h4v)y  '  VJox 

.  „  33^  .  .  38  .  2  hr  3T  _  A 
VDy  w  '  4fcp  3P  -  h8  3  ACi,p  5?  -  0 


(  99) 


The  partial  derivatives  with  respect  to  p  can  be  evaluated  as  follows. 


2  2 
aa  _  3  _  X  =  i_ 

ap  ap  p  o  p 


(100) 


3T  „  3 _  >  p  \  _  1 _ „  T 

ap  ~  3p  pR  ”  p 


35  _ 

ap 


( s-:i ) 


(102) 


When  equations  (100)  through  (102)  are  substituted  into  equation  (.  95), 
the  final  result  is 

2 

h4ux  +  Vy  +  h4  J  !upx  +  vpy)  +  !h2,x.  +  (h3}y  +  u{h4)x  + 

(36) 


For  =  o  >  the  following  result  is  found. 

hlyux  '  <hlyu)x  +  hlyvy  "  ^h1yv^y  +  hiv  + 
+  h2uux  +  hgVUy  +  h3uvx  +  h3vvy  + 

+  VH  -  h4Upx  If-  4  (h4a2»)y  -  h4vpy 

-  h4%ff-hs!AC  •,£*» 


y  ap 


(103) 


The  following  relationships  are  used  to  eliminate  density  derivatives. 


2  2 
3a  _  3  /yPn  _  yp  .  a 

Op  3p  *  o'  2  p 

0 


il  _  3  /£>  _  _£_  _  I 

3c  3p  pR  2  D  p 


t  104) 

{  105) 

(  106) 


Using  the  two  system  momentum  equations  and  equations  (104)  to  (  IC6), 
and  grouping  terms  yields 


-  yuttyx  -  ^(My  -  h2  T  *  h3  ^  4 


+  3  Ch4ux  +  U*h4^x  +  h4vy  +  v(h4)y]  + 

♦  v  4 +  >v*>x  r +  v  ^ +  h4%  r = 

-  A  -2-  (h2(u-upi  +  hj(v-Vp)  +  |cr[h4lT-l) 


-  h. 


(  107) 


Substituting  in  the  Euler  equation  for  z 


2  a  p 

1*1  =  L_  (IR\  =  x 

3X  3X  vp  '  P 


2  a  p 

a«_  =  =  _Jj£.  _ 

3,v  arc'  p 


^  =  p  and  using 


p 


(  108) 


(  1091 


reduces  the  equation  to 

-  yu(h,)x  -  yvOijly  -  n2  h3  ^  + 

+  h«  T  |upx  +  vpy  ‘  aZ(upx  '  v°y>}  - 
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( no) 


-  a2(!i2)x  -  a2(h3)y  =  A  ^  {h2(g-U|))  ♦  ti3(v-vp)}  - 
"  A  (v-1 )  | J  CT[  -1 )  -  h?]j 

iiv:lly,  using  the  system  energy  equation,  equation  (110)  becomes 

IT  PX  +  A  +  ^Vx  +  yv(h1}y  +  a2(h2}x  +  a2(h3}y  = 

=  A  {f  (Y-i )CT[(y-1 )h4  -  h8]  +  yh4B  -  h2(u-up)  -  h3(v-vp}|  (37) 
For  zk  =  up, 

*  h2Aop  ’  h4A?p(B)u  ‘  ^V^x  +  + 


+  Vp^p^x  "  (k6ppup}x  ‘  (h6°pVy  +  ^p  + 
+  VpVx  4  VpVx  =  0 


(  111) 


where 


(B)u  =  -  (y-l)2(u  -  up) 


(  112) 


Rearranging  and  simplifying,  equation  (m  )  becomes 

‘  y°p(h5)x  -  pDUp{h6}x  •  Ppvp(h6)y  +  WVx 
k6^pp^up^x  4  up^cp^x  +  Pp^vp^y  +  vp^pp^yJ  + 

+  h7pp(vp}x  +  h8pp  (hpJx  =  Ap3[h2  ‘  h42(T^)(u“up)  -  h6]  (  1 13) 

Using  the  particle  continuity  equation,  the  final  result  becomes 
h6(up}x  +  h7{Vx  +  h8fhp)x  -  y(h5}v  ’  uP(h6}x  “  vP(h5>y  = 

=  A  |h2  -  h42(Y-l)(u-up)  -  hgj  -  hg  ^  (38 
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F0r  Zk  =  Vp  * 

'  h3App  *  h4App(B)vp  '  +  h5y{ppJy  +  Vp  + 

+  h6(upV  -  (VpVx  "  (h7ppvp)y  +  h7pp{vp}y  +  h7App  + 

-  •.’pVy-0  (114) 

where 

(BL  =  -  (y**1  )2( v-v  )  (115) 

P  H 


Simplifying,  rearranging,  and  substituting  the  particle  continuity  equa¬ 
tion  in  the  same  manner  as  for  z.  =  u  yields 

K  p 


h6(up)y  +  h7(vp!y  {  h8(hp>y  ’  A'y  '  W*  *  vp(h7>y 
=  A  (hj  -  h42(y-l i(v-Vp)  -  h?|  -  h? 

For  z.  =  h^  , 

K  y 

-  b4APp(B)h^  -  (h8Ppup)x  -  (h8Ppvp)y  +  h8  |  ACPp(Tp)hp  =  0 


where 


(Tp-h  *  SE~  ^  *  F* 
y  p  p  c  c 

(B)hp=^[(y-1)fcTp]MY-l)|^ 


(39) 

(  116) 

(  117) 
(  118) 


Substituting  in  equations  (117  )  and  (  118)  and  expanding  terms  yields 


Ppup<h8>x  +  Vp(Vy  +  t,8[pp(up)x  +  “p^p'x  +  Pp(vp>y  +  Wy3 

+  h4flPp<Y-D  fr’-  hjJ»pf  =  °  !  119> 

0  c 
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Substituting  the  particle  continuity  equation  and  simplifying  results 


up'h8^x  f  Vh8-y  =  3? 


|h8  -  (y-l)h4|  +  h. 


!e 

8  y 


in 

(40) 


F°r  zR  =  Pp  , 

hj,A(u-Up)  +  h3A(v-vp)  -  h4AB  +  h5y(up)x  - 

•  (h5yup}x  +  h5y(vp}y  "  fVVy  +  h5vp  =  0  {  120) 


which  reduces  to 


APPENDIX  I 


CHARACTERISTIC  AND  COMPATIBILITY  EQUATIONS 

1 .  GENERAL 

The  method  of  characteristics  is  frequently  used  to  simplify  the 
numeric  calculations  in  supersonic  flow  fields,  where  the  governing 
differential  equations  are  hyperbolic.  In  this  study,  the  governing 
differential  equations  are  a  system  of  hyperbolic,  ouasi-1 inear,  non-homo- 
geneous,  partial  differential  equations  of  the  first  order  as  functions 
of  two  variables.  A  quasi-1 inear  partial  differential  equation  of  the 
first  order  is  defined  as  one  which  is  non-linear  in  the  dependent  vari¬ 
ables  but  linear  in  the  first  partial  derivatives  of  the  dependent  vari¬ 
ables. 

The  method  of  characteristics  develops  compatibility  equations 
which  express  the  dependent  variables  in  terms  of  total  derivatives 
rather  than  partial  derivatives  and  the  characteristic  directions  along 
which  the  compatibility  equations  apply.  This  systen  of  equations  is 
equivalent  to  the  original  set  of  equations  and  is  convenient  to  use  in 
a  finite  difference  form.  The  characteristic  and  compatibility  equations 
can  be  developed  by  use  of  determinants.  However,  the  solution  of  the 
determinants  is  generally  complex  and  an  alternate  method  of  obtaining  the 
characterise  and  compatibility  equations  is  employed  here. 

The  sixteen  equations  which  are  treated  by  the  method  of  character¬ 
istics  consist  of  the  eight  basic  flow  equations  and  the  eight  Euler 
equations.  These  equations  are  rewritten  here  for  convenience. 
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L,  =  Pix  *  pvy  +  upx  +  -'Py  ♦  ^  "  0 

(121) 

4  =  ^x  +  pvuv  +  px  +  A(VU  '  V  =  0 

(122) 

L3  =  puvx  +  pvvy  +  py  +  App  (v  -  vp)  =  o 

(123) 

L4  =  upx  +  vpy  -  a  upx  •  a  vPy  -  ABpp  =  0 

(124) 

P  v 

4  *  Pp<V* +  pp(Vy  +  uP!°p)x  +  Jpi0p>y  +  y  ° 

(125) 

4  ■  Pp“p<“p'x  *  PpWy  •  AoP(u  '  V  =  ° 

(126) 

4  ■  »p“p(Vx "  ppWy '  A°p(v  *  v  =  ° 

(127) 

4  •  yP<Vx  +  ppvp(hp)y  "  f  AC!T  "  V  =  ° 

(128) 

4  =  -P2“x  -  h3vx  •  h4?Px  +  h4  r  Px  +  y(h1,x  + 

+  u(h2)x  +  v{h2)y  -  4  "  0 

(129) 

L10  *  Vy  -  h3vy  -  h4  K  +  P  Py  +  + 

+  u^44 +  v^44  ■  4  =  0 

(130) 

2  2 

.  _  u  j  ^  u  +  h.^"  p  h4  ~~  p„  (44  + 

L11  “  4  x  4  y  4px 

+  (h3)y  +  u{h4)x  +  v(h4)y  -  4  ~  0 

(131) 
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K4  =  A  jj  (Y-l)  CT  [h4(y-l)  -  hg]  -  7^48  - 


“  h2(uV  ’  h3(v‘vp)}  040 

=  *4  |^2  "  h42  ('fl  Hu'up)  "  ht|  “  h6  ^  (141) 

K6  =  fl  {h3  -  V  (Y-Wv-Vp)  -  h7}  -  h7  /  (142) 

K7  '!A§c  (h8-  "4  <1-1 )}  ♦  043) 

Ks  =  A  |h2(u-up;  +  n3(v-vp)  -  h^l  (144) 


It  would  be  desireable  to  separate  the  development  of  method  of  char¬ 
acteristic  equations  into  the  flow  field  analysis  portion  involving 
through  Lg  and  the  optimization  portion  involving  Lg  through  L-jg  since  the 
flow  field  analysis  does  not  depend  on  whether  optimization  is  performed. 
However,  the  presence  of  partial  derivatives  of  the  flow  properties  in  Lg 
through  L-jg  dictate  against  this  approach.  In  reviewing  the  equations,  it 
is  noted,  however,  that  L-|  through  and  Lg  through  L-^  do  net  contain 
partial  derivatives  of  particle  properties  or  multipliers  hg  through  hg. 
Also,  Lg  through  Lg  and  L-j3  through  L-jg  do  not  contain  partial  derivatives 
of  gas  properties  or  multipliers  h-j  through  h^.  Thus  the  method  of 
characteristics  application  may  be  simplified  by  dividing  the  problem  into 
two  segments,  the  gas  property  and  associated  multiplier  equations  and  the 
particle  property  and  associated  multiplier  enuations. 


'V 
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2.  GAS  PROPERTY  AND  ASSOCIATED  MULTIPLIER  EQUATIONS 

A  differential  operator  is  defined  using  arbitrary  functions 
designated  by  o's  as  multipliers  as  follows: 

L  =  o-jL-j  +  O2L2  +  O3L2  +  04  L^  +  cgLg  + 


+  °10L10  +  °nLll  +  c12L12  ::  0 

(145) 

By  grouping  partial  differential  terms,  this  equation  can  be  rewritten 

L  * ft  k  ♦  i  «,}  +  c  | 

*■{*** 

py}  + 

+  G  {"x  +  5  Py}  +  1 

{<».>*  *f  <My}  + 

+  K  {<hZ>x  +  1  <h2>y) 

\  *K  {<h3>x+l  ^3>y} 

+ 

+  p  |(t.4)x*2  (h4)yj 

■  +  R  =  0 

046) 

where  the  coefficients  are 

A  =  po-j  +  p uo2  -  h2  Gg  +  h^o-ji 

(147) 

B  -  pvc?2  *  hgC-jo 

(148) 

C  =  «J03  -  h3cg 

(149) 

D  =  po-j  +  pvg3  -  h3o10  +  h^a-j-j 

(150 

E  =  °2  +  uc4  "  ^4  ^  °g  +  h2  0  G1 2 

(151) 

F  °3  +  V04  '  h4  p  C10  +  h3  0  c12 

052) 

2  2 

6  =  uo1  -  a  uo^j  +  h4  <7g  +  h4  ~- 

cn 

(153) 

H  =  »°1  -  aS  +  "4  4-  =10  +  h4  1 

— 0,, 

3  11 

(154) 
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(155) 


1  =  yoq  +  yuc-jp 

J=  yai0  +  yvo1?  (156) 

K  =  uog  +  o-j-j  +  (157) 

L  =  Vo  g  (158) 

M  *  uo10  (159) 

N  =  Vo1Q  +  g-j-j  +  a^c-|2  (I  GO) 

p=  u0ll  (161) 

Q  =  v0ll  (162) 

R  ■  ‘y'l*  %  (“-up)  °2  +  Acp  !v-vp!  °3  -  ABpp°4  - 

^lcg  "  ^2°10  ~  ^3° 11  ~  ^4°12  (163) 


If  each  of  the  ratios  multiplying  the  partial  derivatives  with  respect 
to  y  are  equal  to  a  value  X  =  dy/dx,  the  differential  operator  reduces  to 

L  =  Adu  +  Cdv  +  Edp  +  Gdp  +  Idh-j  + 

+  Kdh2  +  Mdh3  +  Pdh4  +  Rdx  =  0  (164) 

which  is  the  desired  form  and  is  called  the  general  compatibility  equation. 
Making  the  ratios  each  equal  to  X  results  in  the  series  of  equations 


AX  -  8  =  0 

065) 

CX  -  D  =  0 

(166) 

EX  -  F  =  0 

(167) 

6X  -  H  =  0 

(168) 

IX  -  J  =  0 

069) 

66 


*rr  *  )  j  i 


-W-lAl»ls»jBPS 


3W^4^JWIHP!E" 


-  1  «  n 

«'A  '  U  ‘  V 

/I7n\ 
\  * #  v; 

MX  -  N  =  0 

(171) 

PX  -  Q  =  G 

(172) 

These  equations  are  rearranged  to  group  the  arbitrary  multiplier  terms 
as  the  unknowns. 


pXo^  +  p  (uX  -  v)  a? 

*  ^2a09  ^2°10  ^  ]  =  0 

073) 

-pa-j  +  p(uX  -  v)  o3 

-  h3Xog  +  h3o10  -  h4o^  =  ^ 

(174) 

X02  -  Oj  +  (uX  -  v)  o4  -  h4  1  Xa9  *  h4  i  o10  + 


+  i  (h2X  -  h3)  au  =  0  (175) 

2  a2 

(uX  -  v)  o-j  -  a  (uX  -  v)  a4  +  h4  —  Xcg  - 


- 

>4 

r°io  +  h<s 

j-  (uX  -  v)  0^  =  0 

(176) 

yAOg 

-  ya10  +  y(ux  • 

-v)  oK  =  0 

(177) 

(uX 

-v) 

og  +  Xc11  + 

2 

a  Xc^2  =  0 

(178) 

(uX  - 

-V) 

°io  *  an  * 

a2o12  =  0 

(179) 

(uX  ■ 

-V) 

°11 =  0 

(180) 

If  this  system  of  equation's  is  to  have  a  solution  other  than  the  trivial 
solution  with  the  arbitrary  multipliers  equal  tc  zero,  the  determinant 
cf  the  coefficients  must  equal  zero.  The  determinant  is  shown  in  Figure 
13.  Expanding  the  determinant  results  in 
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\UA  - 


..\4 

v; 


2 

■»  i« 


v  y/\ 


Omi\ 
kVI  f#\ 


n  8i  i 


This  equation  is  now  solved  for  X  tc  give  the  characteristic  equations 
along  which  the  compatibility  equations  will  be  valid. 


2  /? 

^7 


- 1 


(42) 

082) 


Equation  (42)  is  repeated  four  times  in  the  solution  of  the  determinant, 
indicating  that  four  compatibility  equations  will  apply  along  this  char¬ 
acteristic  direction.  Each  of  the  equations  represented  by  equation  (182) 
is  repeated  twice,  indicating  we  can  expect  two  compatibility  equations 
for  each  of  these  characteristic  directions.  For  convenience,  equation 
(182)  is  expressed  in  terms  of  the  flow  angle,  e,  and  the  mach  angle,  a. 

X  =  ^  =  tan(8±a)  (4S) 


Now  that  the  characteristic  equations  are  known,  the  compatibility 
equations  are  determined  from  the  general  compatibility  equation,  equa¬ 
tion  (164).  First,  the  arbitrary  multipliers  in  the  coefficients  must  be 
determined  or  otherwise  eliminated.  This  is  accomplished  by  substituting 
each  of  the  characteristic  equations  into  equations  (165)  through  (172). 
When  X  -  v/u  is  used,  the  following  series  of  equations  result: 


0(h  "  h2  ^9  +  h2°10  +  h4  i  °11  =  C  083) 

-pcj  *  h3  +  ^3°ig  _  h4cn  =  0  084} 

u°2  '  °3  ■  h4  o  S°9  *  h4  ^10  +  o  (h2  I  •  V'k’0  <185> 
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&V2?»5iS 


-2  _2 
<L  v  „  .a 


h4  “p-  i  °9  - 

h4  %  °10  =  0 

(106) 

Z°s  -  °10  = 

0 

(187) 

Jn  +  °i2 

=  0 

(183) 

°11  +  a*  °12 

=  0 

(189) 

0  =  0 

(190) 

Only  four  independent  equations  result  from  this  set. 
-h 

°1  "  h4  p  °12 

°3  =  Xo2  +  p  (h2X  '  V  °12 

°10  =  Xo9 
2 

a-j-j  a  c^2 


(191) 

(192) 

(193) 

(194) 


Four  a's  must  be  taken  as  arbitrary.  Since  does  not  appear  in  the 
four  independent  equations,  it  is  selected  as  one  of  the  arbitrary  a's. 
The  others  selected  are  Og,  ag  and  c-^. 

Now  the  general  compatibility  is  rewritten,  grouping  terms  contain¬ 
ing  each  of  the  four  arbitrary  a's  together  to  yield 

■jpudu  +  ovdv  +  dp  +  Ap„  (u-u Jdx  +  Ap„  (v-v  )dyl  c0  + 

l  p  p  p  p'  2 

+  judp  -  a^udp  -  ABppdxj  + 

+  |»hydu  -  h3dv  -  h^,  ^  dp  +  — -  dp  +  ydh^  + 
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+  udh2  +  vdh^  -  K^dx  -  Kgdyj  o^  + 

'  1  *2 
+  j(vh2  -  uh3)dv  +  h2  j-j-  dp  -  uh4  —  (y-1)  dp  + 

+  yudh-j  -  a2udh4  +  (h^a2  y  +  a2K3  -  K^)dx  + 

+  A  (v-vp)  (h2dy  -  h3dx)j  o12  =  0 


095) 


Since  the  o's  are  arbitrary,  their  coefficients  must  each  equal  zero, 
resulting  in  four  compatibility  equations  which  apply  along  the  charac¬ 
teristic  direction  A  =  v/u, 

pudu  +  pvdv  +  dp  =  -App  j(u-Up)dx  +  (v-Vp)dyj  (1£ 

2 

udp  -  a  udp  *  ABppdx  (IS 


-  h2du  -  h3dv  -  ^4  1  dp  +  h4  —  dp  +  ydh^  + 


+  udh2  +  vdh3  -  fCjdx  -  «2dy  =  0 
1  a2 

(vh2  -  uh3)dv  +  h2  ^  dp  -  uh4  —  (y-l)dp  + 

+  yudh^  -  a2udh4  +  (h^a2  y  +  a2K3  -  K^)dx 


096) 


+  A  ££  (v>vp)  (h2dy  -  h3dx)  =  0 


Along  the  characteristic  directions  given  by  X 
lowing  equation  is  valid. 


097) 


X2(u2  -  a2)  -  2uvA  +  (v2  -  a^)  =  0 


-  tan  (0xa),  the  fol- 


098) 


The  quantity  uX  -  v  cannot  be  zero,  therefore  must  be  according  to 
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from  equations  (178)  and 


(179),  Og  and  o1(0  are  eliminated  from  the  remaining  equations  and  o-j 
determined  from  equation  (176).  This  is  used  to  eliminate  o-j  from 
equations  (173)  and  (174).  Finally,  using  equation  (198)  to  reduce  the 
results,  the  eight  equations  are  reduced  to 


Xa^  _  1  f,  .  Xa^  u  1 

°2  TuX^vT  °4  ~  p  \h2  +  TuX^vJ  h4J 

°12 

(109) 

2  .  ,  2 

°3  (uX-v)  a4  1  p  {(uX-v)  h4  4} 

c12 

(200) 

'°2X  +  °3  *  °4  (uX-v)  +  1  l(uX-v)  h4  + 

hgX  -  h. 

3}  a12 

(201) 

°1  =  a2  °4  +  l  V2  °12 

(202) 

a12  {  0  }  =» 

(203) 

a  -  a2  A 

°9  (uX-v)  °12 

(204) 

a2 

a10  (ux-vi  a12 

(205) 

°11 ■  0 

(206) 

Equation  (201)  may  be  obtained  by  multiplying  equation  (199)  by  -X  and 
adding  to  equation  (200).  Thus,  two  of  the  above  equations,  (201)  and 
(203)  are  not  independent.  Therefore,  two  of  the  o's  must  be  arbitrary. 

and  are  selected.  Substituting  the  six  independent  equations  in- 
to  the  general  compatibility  equation  and  setting  the  coefficients  of 
and  to  zero  results  in  the  following  two  equations. 


* 
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{<»2-SH  du  + 


/as!i  c>v  + 1-  4a< 

(uX-v/  ' 


l  uA-v  tu}dP 

t  {aus.  .  A;^  (U-UD)  +  Apn  (v-vn)  - 


p  UX-V 


-  ABp, 


i  dx  =  0 


PJ 


jfyi2  -  Uh0  -  rr^rr  h„  +  ~~  h„  du  + 


(207) 


2  ux-v  "4  T  ux-v  "2 
2  .2 


v  "4  ’  u,,3  T  HFv  "3}  dv  +{f“  (uh< 


h-  '  «h,  +  ^Ar 
2,  %  (  ..,2. 


a2x 


4  ux-v  “4 
2 


} 


dp  + 


fu  a 
\ux- 

{-y  ^  ♦  «y}  ik,  -  { S0  - a2}  dh2 +  (spr)  dh3 

{7  V2  - A  p£  (h2 +  h4>  (u-uP> + 


+  *^<3P7  h4-A3)  <»-V  +  Ki^- 


-  K, 


,1 


K.  >  dx  =  0 


'2  ux-v  4  j 

After  some  manipulation  including  the  use  of  equation  (198)  and  the 
relationships 

u-v  -  lii-  -  a2>  =  ±  cot  a 


(208) 


(209) 


u  X  -  uv  -  a  x 


■g—  =  +  tan  a 


(210) 


the  compatibility  equations  become 

0  «2  2 

a“  [vdu  -  udv]  ±  —  cotctdp  =  — ~  (udy  -  vdx)  - 
P  » 


-  A  ^-|B(udy  -  vdx)  +  a2  [(u-up)dy  -  (v~vp)dx]l  (21) 
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12- 

+  h^dv  +  -  h4  (dp  -  a  dp)  -  ydh,  +  tana  (vdhg  -  udh^)  = 


±tana  jh3  y  dx  +  A  ^  [h3  -  h42(Y-l)(v-vp)  -  h7]  dx  - 


h2  y  dy  ‘  A  ^  [h2  -  h42(y-l )  (u-up)  -hg]  dy  + 


+  (udy  -  vdx)  A  (Y-l)  [|-  CT  j(Y-l)  h4  -  hgl  +  h4B]j 

(46) 

These  equations  represent  the  four  equations  which  apply  along  the  Mach 
lines. 

3.  PARTICLE  PROPERTY  AND  ASSOCIATED  MULTIPLIER  EQUATIONS 

The  differential  operator  in  the  case  of  particle  properties  and 
associated  multipliers  is 


L  =  °5L5  +  a6L6  +  °7L7  +  °8L8  +  °13l13  +  °14L14  +  °i5L15  +  °16L16 

(211) 

or,  regrouping  by  partial  differential  terms, 

L  '  A  K>*  *  !  '“ply}  *  C  {'Vx  +  !  <Vy}  + 

+  E  {%K  *  I  'hp)y}  +  G  {Vx  +  I  'pply}  + 

+  1  (<hS>x  +  f  (h5>y}  +  {'Vx  +  I  <h6>y}  + 

+  M  {'Vx  +  I  'Vy}  +  P  {'h8>x  +  P  'Vy}  +  R  =  0  '212> 
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where  the  coefficients  are 


A  =  Ppa5  4Ppup°6  +  h6°13 

(213) 

B  =  Ppvp°6  +  h6°14 

(214) 

C  =  Ppup°7  +  h7a13 

(215) 

D  =  ppOg  +Ppvpa7  +  ^7°14 

(216) 

E  =  °pup°8  +  h8°13 

(217) 

F  *  Ppvp°8  +  h8°14 

(218) 

6  =  up°5 

(219) 

H  =  vp°5 

(220) 

I  =  -y  o13  +  upo16 

(221) 

J  =  '^14  +  vp°16 

(222) 

K  =  -upa13 

(223) 

L  =  "vp°13 

(224) 

M=  -upo14 

(225) 

M  =  -V  n, . 

P~  l1* 

(226) 

p  =  up°15 

(227) 

<•'  vp°15 

(228) 

R  *  °5  -  «Pp(u*“p)  °6  ‘  %  (v'V 

°7  "  f  ^_Tp^  °8 
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(229) 


-  K5°13  '  K6°14  -  K7c1S  *  K8°16 
The  general  compatibility  equation  is 

L  =  Adu  +  Cdv  +  Edh  +  Ghp  +  Idhc  +  Kdh,-  +  Mdh7  + 
p  p  p  p  5  6  7 

Pdiig  +  Rdx  =  0  (230) 

Setting  the  ratios  which  multiply  the  partial  derivatives  with  respect 
to  y  equal  to  x,  and  grouping  terms  by  the  o's  results  in 


°pAo5  +  W-V  °6  +  h6X013  "  ,l6°14  =  0 

(231) 

■Op05  J  W'V  °7  +  h7Xo13‘  h7°14  =  0 

(232) 

Pp(upx"vp)  °8  +  h8Xo13  *  h8°14  =  0 

(233) 

('JpJ.-vp)  a5  =  0 

(234) 

-yXo,3  *  yau  +  (upX-vp)  o,6  =  0 

(235) 

‘  V*V  °13  =  0 

(236) 

-  <y-vp)  0,4  -  0 

(237) 

VV  o,5  -  0 

(238) 

The  determinant  of  the  coefficients  is  presented  in  Figure  14.  Expansion 
of  the  determinant  results  in 

(upX-vp)8  =  0  (239) 
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FIGURE  14  DETERMINANT  FOR  FINDING  PARTICLE  CHARACTERISTIC  EQUATIONS 


(22) 


Thus ,  tnc  characteristic  caustic"  is 


>  =  JL 

‘  dx  up 


and  eight  compatibility  equations  should  apply  along  that  direction. 
Substituting  the  characteristic  equation  into  equations  (231)  through 
(238)  results  in  the  following  equations. 


pptf  a5  +  h6  u^  °13  "  h6°14  -  0 
P  P 


(240) 


~pp°5  +  h7  °13  “  h7°14  =  0 


(241) 


h8  ^  °13  -  h8al 4  "  0 


(242) 


a5  (0)  =  0 


(243) 


°13  '  °14  *  0 


(244) 


a,,  (0)  =  0 


(245) 


a™  (0)  =  0 


(246) 


o15  (0)  =  0 


(247) 


which  reduce  to  only  two  independent  equations. 


°5  =  0 

°14  =  °13 


(248) 


(249) 


55^ 


VT 1 ' 


t^J,.  .JJ'JUJJ-  IIIAUil  ■W?I^J.-!4|i.i:ppi 


Setting  the  coefficients  cf  the  arbitrary  multi  pliers* 


&■  r 


u' 


and  ^g,  in  eouation  (230),  the  general  compatibilitv  «auatior  £rf 
particle  properties  and  associated  multipliers*  equal  to  zero  results 
in  the  following  compatibility  equations  along  the  particle  streamline: 


upduP  =  A  (u~up)  dx  '23) 

updvp  =  A  (v~vp)  dx  (24,‘ 

updhp  -  I  «  <T  -  y  d*  <Z5> 

updh5  =  A^h2  ^u_up^  +  h3  ^v‘vp^  “  h48^  dx  (48> 


h,du„  +  h7dv  +  hQdh„  -  ydhc  -  u  dhc  -  v  dlx,  - 
op  7p  8  p  5  p  6  p7 

=  A  jh2  -  h42{y-l)  (u-up)  -  hg|  dx  -  hg  dx  + 

♦A  -jhj  -  h42(Y-’)  (»-»p)  -  l),|  dy  -  h?  dy  (49) 

w  J  f  C"8-  (Y-1)  h4]dx  +  h8^  a*  m 

r  c 

A  deficiency  of  two  compatibility  equations  exists.  None  of  the  equations 
have  derivatives  of  the  particle  density  function,  Pp.  This  deficiency 
is  corrected  by  the  use  of  a  numerical  integration  scheme  to  determine 
the  pp.  The  other  deficiency  prevents  the  calculation  of  the  multipliers 
hg  and  hy  since  derivatives  of  these  variables  appear  only  in  equation 
(49).  This  deficiency  is  corrected  by  using  a  numerical  evaluation 
scheme  to  evaluate  two  of  the  Euler  equations,  (38)  and  (39),  which  per¬ 
mit-  writing  them  in  a  compatibility  like  form  for  the  evaluation  of  hg 
and  hy  along  particle  streamlines.  Equation  (49)  is  not  used,  being 
replaced  by  the  two  compatibility  like  equations. 
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APPENDIX  III 
THE  CORNER  CONDITION 

When  the  end  points  of  an  extremum  are  allcv/ed  to  vary,  the  varia¬ 
tion  must  result  in  no  change  in  the  function,  I,  if  the  functional  is 
at  either  a  maximum  or  a  minimum.  In  the  calculus  of  variations,  this 
restriction  is  expressed  by  the  corner  condition  as  given  by  Mi el e  (17). 

3H,  •»  /■3H,  \ 

A  ]H1  -  r  j’Sx  4  A  Usy  =  0  (5T) 

where  A  denotes  the  difference  between  the  quantity  in  braces  evaluated 
on  each  side  of  the  comer  point  and  6x  and  6y  signify  small  variations 
in  x  and  y.  is  the  integrand  of  the  line  integral  portion  of  the 
functional .  I . 

H1  =  (p  -  P0)nn'  +  +  C2no(un'  -  v}  (31) 

At  corner  A,  the  corner  ooint  is  fixed  and  the  variables  are  fixed 
because  the  nozzle  contour  is  fixed  upstream  of  ooint  A.  Thus,  all 
elements  cf  the  corner  condition  are  zero  and  the  comer  condition  is 
satisfied. 

H1  is  not  apnlicable  at  point  B  since  H-  has  values  along  the 
nozzle  contour  only.  Therefore,  the  value  of  is  identically  zero  on 
both  sides  of  corns**  B  and  the  corner  condition  is  satisfied  identically. 

At  po^nt  C,  the  corner  point  is  not  fixed,  so  5x  and  <Sv  must  remain 
arbitrary.  Therefore  the  coefficients  of  ox  and  5y  must  each  be  zero. 
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MHi 


(250) 


A 


(251) 


Since  H-j  is  exactly  zero  along  the  boundary  8C,  the  quantities  in 
braces  are  identically  zero  on  that  side  of  the  comer.  Therefore,  the 
quantity  in  braces  irast  also  be  zero  at  ooir.t  C  when  evaluated  along  the 
line  AC. 


(p  -  P05nn*  *  C-jG  +  C2*ip(un'  -  *)  - 

*  f  1 

-  y  |(p  -  nQ)n  +  C^.  +  C^ncu^  =  0 


(p  *  PQ}n  +  C^G^i  +  Cgncu  =  0 


(252) 


(252) 


Equation  (252)  is  reduced  by  aoolyinq  equation  (253)  and  the  eouation  of 
the  streamline  at  the  nozzle  contour,  to  obtain 


(p  -  P0)nn'  *  C-jG  *  G 

or,  v'earranging  to  determine  the  value  of  C^,  -nd  using  n!  *  v/u, 

MoVc 


(254) 


‘l  =  " 


(54) 


Substituting  the  value  of  C|  into  eouation  (253)  and  rearranqinq  to 
determine  the  value  of  Cg  yields 

if-cte  6 , .  t  )n 

u,G^  n  'hc  •  o'  c 


V 


c  c 


(255) 


Vcuc 
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or 


(p  -p  )  ■il  -  u  G 
VFc  Fo  \ _ c  c 


c  v 
Fc  c 


(55) 
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APPENDIX  IV 

THE  TRANSVERSALITY  CONDITION 


1 .  GENERAL 


For  an  optimum  nozzle  contour,  the  thrust  must  be  at  a  maximum,  and 
therefore  cannot  change  with  slight  variations  of  the  variables  along 
the  boundary  lines.  Calculus  of  variations  expresses  this  condition  in 
the  transversal ity  equation,  written  here  form  dependent  variables. 


£  <?x  +  n  6y  +  I 
k=! 


(256) 


where  x  and  y  are  the  independent  variables,  zk  (k  =  1  ....  m)  repre¬ 

sents  the  dependent  variables,  and 

''  m 

5  +  y'  [K2  •  kI?k  (K2V  -  *'E  - 

m 

-  [  z'  E  (K.,z.)  (257) 

k=l  K  1  K 

m  m 

n  =  [Kp  -  [q.  (K2)  3  -  y‘  [qk(Ko)n  +  E  (K,,y)  (258) 

£  k=r  ‘  qk  k-r  6  pk  1 


ck  =  ~{h\  +  *  oyp  +  E  (Krzk}  (k  =  1 . m)  (259) 


E  <K;.y)  =  (K,)y  -  57  {(*,),.} 


E  (K,,2k)  *  (K,)z^  -  !;•  {(K,)2,} 


(260) 


(261) 
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(262) 


,«  =  & 


z‘  =  — - 
k  dx 


\  "  3x 


qk  “  3y 


n  +  v  'q 

'  4k 


(263) 


(264) 


(265) 


This  transversal ity  equation  is  based  on  defined  as  the  integrand  of 
the  line  integral  of  the  functional  1,  with  the  line  integral  taken  in 
the  positive  sense  around  the  region  ABC.  K2  is  the  integrand  of  the 
area  integral  of  the  functional  I.  Therefore,  in  the  present  case 


K1  =  “  (p-Pq)™'  ”  C1G  "  Cgnpvun'-v)  along  CA 


(266) 


K1  =  0 


along  AB  &  BC  (267) 


(2  H2  =  Jfih 


The  integrand  K2  is  equal  to  zero  and  K-j  does  not  contain  any  part¬ 
ial  derivatives  of  the  dependent  variables  with  respect  to  x  or  y. 
Therefore  equation  (256)  may  be  rewritten  as 

8  s  % 

y,  »kV  r  e  y,  (Ki)zj 

+  {  l  qkWk  +  E  (Kj  ,y)|  6y  + 


wnere 


+Jll -Mk  +  (Kl>zj  =  ° 


wk =  ■  y’ 


(268) 


(269) 
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Along  the  boundary  line  AB,  all  dependent  and  independent  variables 
are  fixed  since  AB  is  a  right-running  characteristic  from  point  A  and 
the  noz2le  contour  upstream  of  point  A  is  fixed.  Therefore,  the  varia¬ 
tions  <5x,  6y  and  Szk  are  all  zero  and  the  transversal ity  condition  is 
satisfied. 


3.  THE  TRANSVERSAL ITY  CONDITION  ALONG  BC 

Along  the  boundary  line  BC,  is  identically  zero  so  equation  (268) 
reduces  to 


8  8  8 

I  P..W.6X  +  y  q.W.  5y  +  £  W,,6z.  =  0  (270) 

k~1  K  *  k=l  K  K  k=l  K  K 


The  coefficients  of  each  of  the  variations  must  equal  zero  if  the  nozzle 
contour  is  opt1':  urn  since  the  variations  are  arbitrary  along  BC.  There¬ 
fore,  each  of  the  Wk  must  be  equal  to  zero  to  make  the  last  term  zero. 


Wk  =  (K2}qk  ~  y,(K2)Pk 


=  0  (k=l , . . . ,8) 


(271) 


Setting  the  coefficients  of  the  variations  of  gas  properties  to  zero  re¬ 
sults  in  the  following  equations  applicable  along  the  boundary  line  BC. 


n3  -  y'h2  +  h4  (v-y'u)  =  0 

(56) 

Hjy  (v-y'u)  -  h4a2  (v-y'u)  =  0 

(57) 

h2  (v-y'u)  -  h-jyy'  =  0 

(58) 

h3  (v-y'u)  +  ^y  =  0 

(59) 

Variations  in  particle  properties  are  applicable  only  in  the  region 
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where  particles  are  present.  Thus,  setting  the  coefficients  of  the  vari¬ 
ations  of  particle  properties  to  zero  results  in  the  following  equations 
which  apply  along  the  portion  of  the  boundary  between  points  B  and  E. 


-h6pPvp  +  h5»'pp  *  V  Vp  “  0  (272) 

•h5J,pp  -  h7ppvp  +  V,pp“p  °  0  (273) 

-h5yVp  t  h5yy'up  -  o  (274) 

-Vpvp  +  VppuP  ’  0  (275) 

which  reduce  to 

h5(y'up  “  vp)  =  0  (60) 

h6(y‘up  ~  V  +  h5yy*  =  0  (61) 

h7(y,up  -  vp)  -  h5y  =  0  (62) 

h8(y’up  -  vp)  =  0  (63) 


The  variations  of  x  and  y  do  not  provide  any  additional  conditions. 
Both  the  coefficients  of  6x  and  6y  go  to  zero  when  is  set  to  zero. 
Thus,  the  entire  transversal ity  equation  is  satisfied  along  BC  when  equa¬ 
tions  (56)  through  (63)  are  satisfied. 

4.  THE  TRAHSVERSALITY  CONDITION  ALONG  CA 

Along  the  nozzle  contour  CA.  no  particles  are  present,  and  thus  var¬ 
iations  of  partic'.  :•  properties  are  not  applicable.  Using  equation  (268) 
ano  setting  the  coefficients  of  gas  property  variations  equal  to  zero 
results  in  the  following  equations. 


h2(un’-v)  +  h^nn*  -  C2nn*  =  0 

(276) 

n^n  -  h3(un*-v)  -  C2n  -  0 

(277) 

h-jn(un'-v)  -  h4a2(un*-v)  +  C2n(un*-v)  =  0 

(27C) 

h3  "  n<h2  ’  h4(un,~v)  +  ’in*  +  C^Gp  =  0 

(279) 

t' ong  CA,  the  gas  streamline  equation  uq'-v  =  0  applies  and  equations 
(2/6)  through  (279)  reduce  to 


h,  =C2 

(72) 

uh3  -  vh2  +  vq  +  uC-jGp  =  0 

(73) 

Setting  the  coefficient  of  the  variation  6x  in  equation  (268)  to 
zero  yields 

ux  [h2pv  -  h^pn'  -  n2pun']  +  vx  [h-jpn  +  h3pv  -  h^uo']  + 

+  Px  Ch3  -  h2n*  +  h^v  -  h4un’]  +  px  [h^v  -  h^ur,'  - 

-  h4a2v  +  h4a2un‘]  +  (p-PQjn’n*  +  C3 +  C2p(uo‘-v)n'  - 

-  n*  [(p-pQ)n  +  C-,6^.  +  C2nou]  -  [-C^pn1]  - 

-  gj  tc2np]  ♦  §  [nn'  +  C,6p]  -  jjf  [C2n(un'-v)]  =  0  (  280) 


Using  the  definition  of  a  total  derivative. 


dF  s  3F  ,  3F 
dx  '  3x  n  3n 


(281) 
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and  equations  (276)  through  (279)  and  the  equation  of  a  streamline, 
un'-v  =  0,  equation  (280)  is  reduced  to 

dC2  dG  . 

-nou  dl"  +  ClGn  ~  C1  “dx~  +  ClpnGp  “  pxn  " 

-  c2  {Slnfiyl  .  n,npUri  +  npv^j-  0  (282) 

Using  the  continuity  equation 


npv^  -  -  npux  -  nupx  -  nvp^  -  pv  = 


npux  -  nupx  -  nn'up^  -  n'up  (283) 


and 


=  npux  +  nupx  +  n‘npun  +  n'nup^  +  n'pu  (284) 

the  term  in  the  brackets  j  j  goes  to  zero.  Further,  with  no  particles 
present,  the  momentum  equations  can  be  written  as 


Px  =  -  puux  -  pvun  =  -Pb  £ 


"n  =  -  puvx  -  pvvn  =  -pu  37 


(285) 

(286) 


which  together  with  the  relationship  C9  =  h-j  reduce  equation  (282)  to 


dh 


1  _  du  ,1  L  d  #dG 


36 \  dv\ 


dx  dx  npu 


Jr  U  /Ola  N  uv 

\Gn  ■  37  V0  '  pu  3p  37/ 


(74) 


Setting  the  coefficient  of  the  variation  <5y  in  equation  (268)  to 
zero  yields 
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un  [h2pv  -  h^npo *  -  h2pun‘]  +  v^Ch^np  +  h3pv  -  h3pur,']  + 

+  p^[h3  -  h2n‘  +  h^v  -  h^un']  +Dn[h-|nv  -  h^nun'  - 
-  h4a2v  +  h4a2un‘]  -  (p-P0)n*  -  -  C2p(un'-v)  + 

+  3x  C(P“PQ)n  +  ciGr|l  +  C2nou^  =  G  (287) 

Using  the  streamline  equation  and  equations  (276)  through  (279)  reduces 
equation  (287)  to 

dC2  dGn« 

-npu  +  C]Gn  -  C,  -gjp  +  C-jP^Gp  -  pxn  - 

-  C2{^^  -  n*ncun  +  opvn}  =  °  (282) 

which  is  the  same  result  as  obtained  for  the  variation  in  x.  The  trans¬ 
versal  ity  equation  (74)  therefore  accounts  for  variations  in  both  x  and 

y- 
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